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Some Theorems and Procedures for 
Loop-F ree Routing in Directed 
Communication Networks 


By R. MAGNANI 


(Manuscript received November 15, 1967) 


This paper examines two general methods for specifying directed routing 
patterns 1n communication networks. Hierarchical routing, as currently 
used in the toll network, is such a directed routing pattern. However, it 
ts only one member of a large set of possible routing strategies that can be 
realized by storing, in each office, a list of outgoing trunk groups and an 
order-of-chotce for these groups for each received call address. The general 
class of routing strategies is defined by this method of realization, subject 
to the constraint that routing patterns be loop-free. The paper discusses 
procedures for generating loop-free patterns, for detecting whether or not 
a gwen pattern ts loop-free, and for specifying ‘‘good” patterns from the 
large number which are realizable. 


I, INTRODUCTION 


The fundamental problem which besets people concerned with the 
design of communication networks is how to provide a network which 
is, at once: 


(t) Of sufficient routing capability to allow any two users to be 
connected with a high probability of success. 

(2) Economical in its use of transmission facilities and switching 
centers. 
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(wit) Capable of surviving extensive natural or man-made damage. 

(iv) Adaptable to changing traffic patterns and overload situations. 

(v) Capable of being engineered and implemented in small sections, 
over a period of years, by many different people. 


This is a problem of such complexity that, with the present state of 
the theory, it must be attacked piecemeal. 

This paper deals with a small but important section of the problem. 
It considers some of the topological properties of communication net- 
works and examines a class of alternate routing strategies from a gen- 
eral point of view. Our purpose is to state rules which will allow the 
orderly production of routing patterns, for arbitrary networks, by a 
computer. We approach the problem in three stages: 


(t) Several simple rules are stated for producing “loop-free” rout- 
ing patterns. 

(72) The rules are “generalized” to allow the proof of some theorems 
about the extent of their application. 

(ww) A more limited but practical statement of the rules is presented 
followed by several heuristic procedures, based on these rules, which 
can be used to specify “good” routing patterns from the large number 
which can be generated. 


II. BACKGROUND 


What is fundamental to the routing process as it is practiced today 
in the telephone network? Certainly one of the answers is that each 
office stores a list of outgoing trunks and an order of choice for those 
trunks for each possible call address that can be received. We can 
think of the aggregate of these lists as constituting a “routing map” 
which has been distributed among many offices. The “map” which is 
currently stored in the telephone network realizes the hierarchical 
routing plan. 

But suppose we wish to examine strategies which do not require a 
“hierarchy” of offices? Is there a way to realize a general class of 
routing maps which will allow offices to be of equal rank and which 
can be implemented in the same fashion as the current hierarchical 
plan? The answer is that such a class of routing maps does exist and 
that, indeed, the present hierarchical map is simply one member of 
the set. To see this, consider the simple network of four offices shown 
in Fig. 1. This may be an entire network or some subset of a much 
larger network with the connecting trunk groups omitted. For the 
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present, we will assume the trunk groups shown are all two-way (that 
is, contain trunks that can be seized from either end) and that rout- 
ing between the offices is subject to the following constraints: 


(¢) No routing control information is passed between offices. 

(iz) Offices do not check for shuttle.* 

The resulting network is a fair approximation to a subset of the 
present day commercial network.+ Routing between offices is accom- 
plished as follows : 


(1) Each office is assigned a unique address such as NNX or NPA-NNX. 

(ii) Upon receiving a call request, an office checks to see if it is the 
destination office. If it is, the call is counted as a success (although 
in practice the called subscriber’s line must still be checked). If not, 
the office consults a routing table and, on the basis of the received 
NNX, selects an outgoing trunk group. 


Fig. 1— A four-office network. 


(iit) If a trunk in the group is available, it is seized and the called 
address is passed over it to the next office. At this point we return to 
step (22). 

(iv) If no trunks are available, the office again consults its routing 
table, to find an alternate trunk group, and returns to step 2. The 
process continues until all alternate trunk groups have been tried. 

(v) If no trunks are available in any of the alternate trunk groups, 
the call is blocked, and the caller is so notified. 


For a particular network and destination office, this process may be 
conveniently summarized on the graph which represents the network; 
for routing to a particular office, we assign each trunk group in the 
network a direction and an order of choice out of the office in which 
it originates. For example, if office 4 were the destination office in the 


* Shuttle refers to routing a call out over the same trunk group on which the 
call arrived. : 
+ With one-way trunk groups omitted. 
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network of Fig. 1, we might summarize routing to this office by the 
use of Fig. 2. 

Here we are to understand that calls must route in the directions 
shown and that trunk groups are to be selected in the given order 
(beginning with 1). The routes between office 1 and office 4 can easily 
be seen to be: 1-2-4, 1-8-4, and 1-38-2-4, where the numbers rep- 
resent office numbers and the routes are listed in the order they will 
occur. Similarly, offices 2 and 3 can be seen to have routes 3-4, 3-2-4, 
and 2-4, 

Fig. 2 represents routing from all offices to office 4 and will be said 
to be a directed routing pattern to office 4 on the network of Fig. 1. 
When such a routing strategy is followed, the way in which a call 
routes from a particular office is independent of the past history of the 
call; this is characteristic of routing in the present DDD network. To 


DESTINATION 
OFFICE 





Fig. 2— A directed routing pattern to office 4. 


completely specify routing in the network of Fig. 1, four directed 
routing patterns are necessary, one with each office as destination. 
The direction and order of choice assigned to the trunk groups will 
vary from pattern to pattern depending on the destination office. To 
see that the hierarchical pattern in use today is of this type, we need 
simply draw the pattern. (See Fig. 3.) 


III. ROUTING PATTERNS 


Let us examine some simple rules for constructing such patterns 
and adopt the standard graph theory terminology: “branch” or “link” 
for trunk group, and “‘node” for office. 

Because we assume offices do not check for shuttle or looping, we 
will require the patterns we generate to be loop-free.* A “loop” for 
our purposes is defined as: a set of branches and nodes (not contain- 
ing the destination node) constitutes a loop if we can select any node 


*It has been suggested by J. H. Weber that_a small probability of looping 
may be acceptable if looping can be detected (see Ref. 1), 
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in the set and, by following the directions of the branches, traverse 
each branch once to form a path which returns to the selected node 
(that is, loops must be “directed’’). It is not clear in the case of a 
large network just how one goes about obtaining a loop-free directed 
routing pattern, particularly if the network is nonplanar. To dem- 
onstrate that a systematic procedure is required, we invite you to try 
to draw a loop-free pattern on the network of Fig. 4a. 

A closely related problem is illustrated by Fig. 4b in which we are 
given a routing pattern and asked to determine whether or not it 
contains a loop. In this case, the single loop that the network does con- 
tain may or may not be obvious to you; however, if the network were 
much larger, say 40 nodes, a systematic procedure again would be 
required. 





Fig. 3— Hierarchical directed routing pattern. 


Consider the following two rules for generating directed routing 
patterns in an arbitrary network. 

Rule 1—Select any node in the network (usually the destination 
node) and label all its branches incoming. 

Rule 2—Now select a node which has at least one outgoing branch 
and label all its remaining free branches incoming.* Repeat this 
rule until all branches in the network have been assigned a direction. 

Fig. 5 illustrates the process for a simple 6-node network. The num- 
bers in the node circles represent the order in which rules 1 and 2 are 
applied, beginning with the heavily circled node that is the destination 
node for this pattern. Notice that the process is finished in four steps, 
leaving the two blank nodes with only outgoing branches. We will call 
nodes of this type (the blank nodes) originating nodes, although it is 
assumed that calls routing to the destination node can originate at any 


* Free branches are those which have not been given a direction. 
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node. Nodes with both incoming and outgoing branches will be termed 
tandem nodes (for example, 2, 3, 4). The remaining node, node 1, will 
be called a terminating node and, in this case, it is the destination 
node for the pattern. Indeed, if the rules remain as stated, there can 
be only one terminating node in any pattern, the destination node. 
Unless stations are being considered which are multiple-homed, this 





Fig. 4—A sample pattern. 


is not a limitation.* However, to deal with a multiple-homed situa- 
tion and to allow the proofs of some general theorems, we will remove 
this restriction by restating rule 1: 


Rule 1’ Select any free node in the network and label all its branches 
incoming.+ This rule may be repeated an arbitrary number of times, 
each application creating a terminating node. 


* A station which can be reached from more than one office is said to be 
multiple-homed. A dual-homed station, for example, can be reached from either 
of two offices—in this case, we would want both offices (if not connected) to be 
terminating nodes in the routing pattern. 

+A free node is one with no directed branches. Each application of this rule, 
of course, creates a “trap” for traffic. 
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Rule 2’ Same as 2. 

This revised set of rules will be referred to as backward production. 

Clearly an analogous process exists in the forward direction and 
will be called forward production. 

Rule 1” Select any free node in the network and label all its branches 
outgoing. This rule may be repeated arbitrarily, each application 
creating an originating node. 

Rule 2” Now select a node which has at least one incoming branch 
and label all its remaining free branches outgoing. Repeat rule 2 until 
all branches in the network have been assigned a direction. 


We show later that backward production is the more useful process 
for generating telephone network routing patterns. 


BACKWARD PRODUCTION 
RULE 1: 


BECOMES 


RULE 2: 
BECOMES 


DESTINATION 
NODE 
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RULE 1 
APPLIED 


RULE 2 
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P 
RULE 2 ATTERN 
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Fig. 5— Example of use of backward production. 
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IV. GENERAL THEOREMS 


4.1 Proof of Lemma and Theorems 
We now prove the following lemma and theorems: 


Lemma 1: Any loop-free pattern must contain at least one originat- 
ing node and at least one terminating node.* 

Proof: This can be proven by exhausting the possibilities. Clearly it 
is not possible to have a network consisting only of originating nodes 
or only of terminating nodes. The remaining possibilities are: 


(t) Only Tandem Nodes—If a is a tandem node, it must have an 
outgoing branch to some other node, Bs. Similarly, B must have an out- 
going branch to some node other than a (or we would have a loop), 
say c. This argument proceeds until all nodes in the network have 
been considered. The last node to be considered must connect to some 
previous node since it, too, is a tandem node. Such a connection would 
create a loop. 

(ii) Originating Nodes and Tandem Nodes—lIf a is a tandem node, 
its outgoing branches must connect to another tandem node, since no 
branches can terminate on an originating node. Therefore, the argu- 
ment presented in 2 can be used here. 

(12) Terminating Nodes and Tandem Nodes—If a is a tandem 
node, its incoming branches must originate at another tandem node, 
say B, since no branches originate at terminating nodes. Since B is 
also a tandem node, it must have incoming branches from some node 
other than a (or we would have a loop), say c. Again, this argument 
proceeds until all nodes in the network have been considered. The last 
node considered must have an incoming branch from some previously 
considered tandem node, thus creating a loop. 

The remaining two possibilities (terminating the originating nodes 
only, and all three node.types), each contain at least one terminating 
and one originating node. 


Theorem 1: Routing patterns generated by backward production are 
loop-free. 


Proof: Rule 1’ tells us we may create terminating nodes arbitrarily 
(we must create at least one) in sequential fashion. Since candidates 
for terminating nodes must have all branches free and these branches 

* As this paper was being prepared, the author learned of work by S. L. 


Hakimi in which Lemma 1 and Theorems 1 and 2 are proved in a more formal 
fashion. (See Ref. 2.) 
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are all made incoming, it is not possible to loop through a terminating 
node, nor can there be branches between terminating nodes. Let A 
and ps be terminating nodes and let c be the first nonterminating 
node, with a branch to a (and possibly to s), to which we apply rule 
2’. Since it is not possible to loop through nodes A and B, we may 
disregard the branches to these nodes as far as loops are concerned. 
Then the only branches which can be members of a loop are the re- 
maining (all-free) branches on c. But rule 2’ tells us to make all these 
branches incoming. Therefore, the only way to loop through node c 
is to loop through node a (or B). Since it is not possible to loop through 
A or B, it is not possible to loop through c. Clearly the same argument 
applies at each stage of the process; the only way to loop through the 
present node is to loop through some previously considered node, 
which is not possible. The process ends when all branches have been 
given a direction. At this point, the originating nodes will be seen to 
have been created by applying rule 2’ to all the nodes to which they 
connect. Since it is not possible to loop through originating nodes, the 
pattern must, indeed, be loop-free. 

By a completely analogous proof, it may be shown that the routing 
patterns generated by forward production are also loop-free. 

Now we have two procedures for generating loop-free patterns. The 
question is: What sort of patterns do they generate? We prove: 


Theorem 2: All loop-free routing patterns can be generated by back- 
ward production. 


Proof: This is proven by induction on Lemma 1. Let No be any arbi- 
trary loop-free routing pattern. Then, by Lemma 1, it must contain 
at least one terminating node. For generality, assume it contains two, 
A and s. We will make these nodes evident, leaving a neo? net- 
work, Ny. See Fig. 6. 

In a blank network (that is, a copy of No without branch assign- 





NETWORK. No BLANK NETWORK 


Fig. 6— Construction of a duplicate of No from a blank network—first stage. 
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ments), the terminating nodes a and B may be generated by the ap- 
plication of rule 1’. 

If we were to remove these nodes and their branches from No, we 
would have left the network N;, which must also be loop free. (If Ni 
is not loop-free, then neither is No.) Since N, did not contain terminat- 
ing nodes (all terminating nodes in No were made evident), we must 
create terminating nodes in the process of discarding nodes a and B 
and their branches. That is, in the set of nodes to which a and B con- 
nect, there must be at least one node which becomes a terminating 
node if its branches to A and B are removed. In addition, if there is 
more than one such node, the nodes cannot be connected to each other. 
(Any such connection would make one of the nodes a nonterminating 
node.) We will call these terminating nodes, generated by discarding 
previous nodes and branches, pseudoterminating. Let us assume there 
are two such nodes, c and D, in the network N, and place them in evi- 
dence. See Fig. 7. In the blank network, the outgoing branches from 
c and p (and all other nodes) to A and B, were generated when we 
applied rule 1’ to nodes A and s. If we were now to apply rule 2’ to 
nodes c and pb, in any order, we would generate the nodes c and p in 
the blank network exactly as they appear in the network No. 

If we now remove nodes c and p and their branches from N,, we 
are left with network N2, which must also be loop-free. No had con- 
tained only tandem and originating nodes (all pseudoterminating 
nodes in Ny were made evident). With the removal of branches to c 
and Dp, we must therefore create at least one pseudoterminating node 
in Ng. Let there be two such nodes, & and F, and make them evident. 
See Fig. 8. In the blank network, the application of rule 1’ to nodes 
A and B, and of rule 2’ to nodes c and pb, assigned all the outgoing 
branches from nodes E and F (and all other nodes) to nodes 4, , c, and 
p. Now the application of rule 2’ to nodes & and F in the blank net- 
work will properly assign all the incoming branches to nodes § and r, 
and these nodes will appear as they do in No. 





NETWORK Ny BLANK NETWORK 


Fig. 7— Second stage. 
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NETWORK Noa BLANK NETWORK 


Fig. 8— Third stage. 


This argument may be applied to smaller and smaller networks, 
each time generating the proper nodes and branch assignments in the 
blank network. Since the network No is assumed to be finite, the proc- 
ess terminates in K steps with some network, Nx, which contains all 
the originating nodes in No, now isolated (all branches will have been 
discarded). At this point in our assignments in the blank network, 
we will find that all branches have been assigned and that we have 
generated network No by applying rules 1’ and 2’. 

Again, an analogous proof will show that all loop-free patterns may 
be generated by forward production. 

It is possible to decide whether or not a given network and rout- 
ing pattern contains a loop by a procedure which is a variant of that 
given in theorem 2. (This procedure relies on lemma 1 for its justifica- 
tion.) 

Identify all originating and terminating nodes and remove them and 
all their branches from the graph. Now repeat this step until either: 


(1) Only branchless nodes remain, or 

(it) No originating or terminating nodes can be found. If the net- 
work can be reduced to branchless nodes, it is loop free. If not, at the 
point where no further reduction is possible, the remaining network 
contains at least one loop. 


4.2 Conclusions from the Theorems 


We can now draw the following conclusions: 


(1) The class of routing patterns defined by rules 1’ and 2’, or by 
rules 1” and 2”, is equal to the class of all loop-free routing patterns. 

(wt) Therefore, any pattern that can be generated by forward pro- 
duction can also be generated by backward production. 

(viz) If a pattern is loop free, there is at least one sequence of nodes, 
to which we can apply backward production, that will generate the 
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pattern. (A similar statement holds for forward production.) The 
sequence that will generate the pattern is discoverable by applying 
the procedure given in theorem 2. 

(iv) If N; represents the number of terminating or pseudoterminat- 
ing nodes made evident in the 7" step, the number of ways to generate 
the pattern No, by backward production, is: 


K 
No. of Ways = [J (N; !) 
t=1 
In general, the number of ways to produce a pattern using forward 
production is unequal to the number of ways using backward produc- 
tion. 


V. SYMMETRIC NETWORKS 


5.1 Additional Conclusions 


The theorems and conclusions of Section 1v apply equally well to 
symmetric networks (that is, fully interconnected). In addition, it 
can be shown that, for symmetric networks, the following is true: 


(t) There is exactly one originating node and exactly one terminat- 
ing node in any pattern. 

(77) There is only one way to produce a given pattern using back- 
ward production. (Similarly for forward production.) 

(zt) If we choose a terminating node and apply backward (or for- 
ward) production, we can generate (N-1) distinct loop-free routing 
patterns to the terminal node. These patterns will all be isomorphic 
(that is, the same with a relabeling of the nodes). Indeed, there is 
really only one “abstract” loop-free pattern in a symmetric network, 
no matter which node is the destination node or what order of choice 
is applied, since all patterns can be shown to be isomorphic. 


5.2 Routing in Symmetric Networks 


We would now like to examine routing in symmetric networks both 
as useful in itself and for a bound on routing in incomplete networks. 
We begin by deriving an expression for the number of K-link routes in 
a symmetric routing pattern from a given node to the destination node. 

Assume we are given an N node symmetric network with the first 
node as the destination. Also assume forward production is applied to 
the network, using rule 1 on node N and rule 2 on the nodes N — 1, 
N — 2,-+--+ 3, and 2, in that order. (Equivalently, we could use backward 
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production, applying rule 1 to node 1 and rule 2 to nodes 2, 3, 4, --- 
N — 2, and N — 1, in that order.) Then we may show the following: 


Theorem 8: In a symmetric network of N nodes, the number of K-link 
routes between the Q** node (2 S Q < N) and the destination node (node 1) 


as exactly given by the binomial coefficient C ( Z _ a 


Proof: Consider node N. It will have one branch to the destination 
node, giving a single one-link path. We may write this as Ge a Z . 
Now, any 2-link path must pass through an intermediate node, of which 
there are N — 2. If we can show that each of the ee = ?) selections 
of an intermediate node generates exactly one 2-link path, the number 
of 2-link paths from node N to the terminal node will be o(% ‘i: a 


Let A be one of the possible intermediate nodes. Since we are using 
forward production to generate the routing pattern, each successive 
node considered must have no branches directed toward previously 
considered nodes and must have exactly one branch directed to each of 
the nodes not yet considered. Since node A is considered sometime after 
node NV, there must exist a branch from N to A. But A is considered 
before the destination node (which is last in the process); therefore, 
there must exist a branch from A to the destination node. Thus, there is 
exactly one 2-link route from node N, through node A, to node 1. This 


argument is valid for any of the C @ i. 


there are exactly an it ?) 2-link routes from node N to node 1. 


?) choices for A. Therefore, 


This argument may be generalized for K-link routes. Each K-link 
route requires K — 1 intermediate nodes between node N and node 


1. There are ao = *) ways to choose a distinct set of K — 1 inter- 


mediate nodes. If we let A, , A., A3, --* Ax-; be a particular choice 
of the K — 1 nodes, then there must exist exactly one ordering of these 
nodes which represents the sequence in which rule 2” was applied. 
If the ordering is A, , Az, A3, --: Ax-;, node A, will have a branch 
to A, , which will have a branch to A; , and so on. Since A, is always 
considered after node N, and node Ax-_, is always considered before 
node 1, there will be exactly one K-link path for this choice of K — 1 
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N — 2 


K — ) K-link 


nodes. Since the choice was arbitrary, there are o( 


routes from node N to node 1. 

Now consider node N — 1. All branches from node N were made 
outgoing and therefore are of no use to later nodes for the purpose of 
routing. If we remove node N and its branches from the graph, we 
are left with an N — 1 node symmetric network. In this reduced graph, 
we may consider node N — 1 as the originating node and node 1 as 
the destination node. The sequence in which rule 2” was applied in this 
graph is the same sequence used in the larger graph. Hence, the argu- 
ment presented above applies to this network with N — 1 nodes 
replacing N nodes. The number of K-link routes is therefore 
of @ K D i atl We may proceed to remove node N — 1 and its 
branches from the graph, and so on, generating successively smaller 
symmetric networks and applying the same arguments at each stage. 
In general, from node Q, the number of K-link routes to the destination 


isC(@—?),@< Qs .Nand1 <K <Q — 1). Asa corollary, it 


can be shown that the total number of K-link routes from all nodes 


in the graph to the destination node is given by C Gi 7 t That is: 


(Me )= Egy) 
where 


(2 = q) is zero for K 2 Q. 
It is possible to summarize routing in symmetric networks by using 

Table I. As an example of the information obtainable from the table, 

consider a 6-node symmetric network to which we have applied back- 

ward production in the order: (rule 1) 1, (rule 2) 2, 3, 4, 5, 6. (See 

Fig. 9.) 

This network will have: 


From node 6: one 1-link, four 2-link, six 3-link, four 4-link, and one 
5-link route to node 1 (the destination). 

From node 5: one 1-link, three 2-link, three 3-link, and one 4-link 
route to node 1. 

And so on, reading node I routes from line I. Reading from N=6, 
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there will be a total of five 1-link routes from all nodes to node 1, 
a total of ten 2-link routes, and so on. 

These numbers represent the maximum numbers of K-link routes 
in an arbitrary (not necessarily symmetric) 6-node network. This 
follows from the fact that we may generate the arbitrary network by 
removing branches (and therefore routes) from the corresponding 
symmetric network. 


VI. HEURISTIC PROCEDURES FOR ARBITRARY NETWORKS 


_ How does one go about choosing a “good” routing pattern to a 
given terminal node in an arbitrary network? We can begin by mak- 
ing the following observation. 





Fig. 9 — Six-node symmetric network. 
; Since a practical pattern defines routes to one node, the destination, 
it is reasonable to require “good” routing patterns to have only one 
terminating node type, the destination node. This means we may 
apply backward production as defined by rules 1 and 2; therefore, 
we cannot expect to generate all loop-free routing patterns (which 
require rule 1’ in place of rule 1). This is not a limitation unless we 
are dealing with multiple-homed stations. We will ignore this latter 
case, although the procedures we discuss can be generalized to deal 
with multiple homing. 
_ Now, what is meant by a “good” routing pattern? One with the 
lowest average blocking from all nodes to the destination node? A 
pattern which minimizes blocking from a selected node to the destina- 
tion? One with the smallest average route length? A pattern with the 
maximum total number of routes?* 

To the author’s knowledge no algorithm exists which will guarantee 

* There is, of course, the larger problem of designing a network which realizes 
all of these and is, at the same time, rugged, economical, and so on, as described 


in the introduction. This is a complex problem; its very statement is difficult and 
has been the subject of intensive study, (See Refs. 3, 4 and 5.) 
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any of these criteria. However, it is possible to approach the last two 
criteria by using a heuristic procedure which will generate patterns 
with large numbers of short routes, and which also has the virtue of 
assigning orders of choice to the branches. 


6.1 Generating Patterns with Many Short Routes 
Consider the following method for applying rules 1 and 2 to an 
arbitrary network: 


(1) Select the destination node as the first node and apply rule 1, 
labeling all the branches incoming. Label the originating ends of each 
of these branches the first choice out of the respective nodes. 

(wz) Now, in the set of nodes to which the destination node con- 
nects (these will be called “predecessors” of the destination node), 
select any node and apply rule 2, labeling its free branches incoming. 
Label the originating ends of these branches first choice out of the 
respective nodes, if possible; or, if a first choice already exists (from 
step 2) label the branch second choice. 

(wz) Continue step 2, choosing nodes only from the predecessors of 
the destination node; each time, label the branches n plus first choice 
out of the node at the originating end, where n choices already exist. 
Continue until all the predecessors of the destination node have had 
rule 2 applied to them. 

(tv) Consider the set of nodes which has outgoing branches to any 
node (or nodes) which are predecessors of the destination node.+ 
These may be thought of as second level predecessors of the destina- 
tion node. Apply rule 2 to these nodes until they have been exhausted 
(or until you are exhausted, whichever comes first), each time label- 
ing branches the n plus first choice out of the node in which they 
originate. . 

(v) Identify the third level predecessors of the destination node, 
and so on. Continue the process until every branch in the graph has a 
direction and order of choice out of the node from which it originates. 


Fig. 10 gives an example of the procedure, which is tedious to de- 
scribe, but easy to perform. 

At this point, we can observe that all the paths from the K® level 
predecessors of the destination node have at least K links. We prove 
the following theorem: 


7 The predecessors of any node (or nodes) can be identified without reference 
to branch directions. In this procedure, a predecessor of node A is any node 
connected to node a by an (as yet) undirected branch. If we are seeking the 
pee of a group of nodes, branches between nodes in the group are 
ignore 
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‘PREDECESSORS FIRST LEVEL SECOND LEVEL 
f™N 


RULE 4 
TO NODE 1; 


RULE 2 
TO NODE 2: 


RULE 2 
TO NODE 3: 


RULE 2 
TO NODE 4: 





Fig. 10 — Deriving order of choice in the backward production process. 


Theorem 4: The given procedure creates the maximum number of 
1-link and 2-link routes. 


Proof: Clearly, there is no way to create more 1-link paths to the des- 
tination node than to label all its branches incoming. Consider the 
first level predecessors of the destination node. Among their free 
branches, they may have branches to each other, and branches from 
second level predecessors. 

Each time we label a branch between first level predecessors we 
create one 2-link path to the destination. This is true regardless of 
the direction given to the branch. Hence, the number of 2-link paths 
created this way is fixed, and is exactly the number of branches be- 
tween first level predecessors. If we now discard the branches between 
first level predecessors and consider the reduced graph, it is clear that 
the way to get the maximum number of 2-link paths is to label every 
free branch, on every first level predecessor, incoming. But this is 
exactly the effect of the given procedure. Branches that do not con- 
nect first level predecessors remain free until rule 2 is applied to the 
node; then they are all labeled incoming. It follows that the total 
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number of 2-link paths created is fixed, and is equal to the number 
of free branches on all first level predecessors after the application of 
rule 1 to the destination node. The order in which rule 2 is applied to 
these nodes has no effect on the number of 2-link routes. 

It might seem that this theorem can be extended to show that the 
procedure produces the maximum number of K-link routes (K > 2), 
subject to the fact that K-1 link, K-2 link, ... , 2-link, and 1-link 
routes have been maximized. Unfortunately, one need go no higher 
than 3-link routes to find a counterexample as shown in Fig. 11. 

The heuristic procedure can be improved by eliminating the arbi- 
trary choosing of nodes in step w and in later steps. That is, having 
identified the N® level predecessors of the destination node, we apply 
rule 2 to these nodes in a particular order. 


6.2 Choosing N‘* Level Predecessors 
We suggest this revised heuristic procedure for choosing among N® 
level predecessors: 


(z) Arrange the graph to show the various level predecessors in 
stages. An example is Fig. 12, where higher and higher level predeces- 
sors are encountered as we progress from left to right. 

(iz) Direct all branches between stages toward the destination node. 


(See Fig. 12.) 
(wi) Now consider the first level predecessors, nodes 2, 3, and 4. 





Fig. 11—Counterexample. (a) Pattern generated by heuristic procedure; 
number of 3-link routes: 2. (b) Pattern with maximum number of 3-link routes; 
number of 3-link routes: 4. 
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For each of these nodes, we compile two figures, the number of routes 
provided from the node to the destination, and the number of nodes 
served by this node. Node B is said to be served by node A if there 
exists at least one directed path from B to a. Node 2, for example, 
serves nodes 5 and 9, while node 7 serves none. We take the difference 
of these two numbers (nodes served minus routes provided) and use 
the resulting number as a measure of the need for additional routes. 
For Fig. 12b these numbers may be tabulated as follows: 


No. Nodes No. Routes 
Node Served Provided Difference 





2 2 1 1 
3 3 1 2 
4 3 i 2 


(iv) Now choose the lowest of the difference numbers and apply 
rule 2 to the corresponding node. In this example, node 2 is the choice 
and we label all its branches incoming. (Presumably, it needs the 
least number of additional routes.) 

(v) Node 2 is now removed from consideration and we may restate 
the table for nodes 3 and 4, adding the routes picked up by the 
branches directed into node 2: 


No. Nodes No. Routes 
Node Served Provided Difference 


3 3 2 1 
4, 3 2 1 


In this case, we have equality and so choose node 8 arbitrarily. Node 
4 is, then, the last node in the process and the result is shown in Fig. 
12c. 
(vt) We now move onc stage to the right and consider second level 
predecessors: 
No. Nodes No. Routes 
Node Served Provided _ Difference 


5 1 1 0 
6 2 2 0 
7 0 4 —4 
8 1 4 —3 


This suggests that node 7 is least in need of additional routes and we 
may label all its branches incoming. Restating the table two more 


LOOP-FREE ROUTING 485 


PREDECESSORS : SECOND 
FIRST LEVEL THIRD 


DESTINATION 
NODE oO 





Fig. 12—-Example of heuristic procedure. 


times, we obtain node 8 next and, finally, node 6. The result 1s shown 
in Fig. 12d. 


No. Nodes No. Routes 
Node Served Provided Difference 


5) 1 1 0 
6 2 6 —4 
8 1 8 —7 
5 1 1 0 
6 2 14 —12 


To obtain an order of choice for the branches, we simply apply the 
heuristic procedure for generating patterns with large numbers of 
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short routes in node order 1, 2, 3, 4, 7, 8, 6, and 5. (See Section 6.1.) 
The result, identical to that in Fig. 12d, is shown in Fig. 18. 

This method yields routing patterns in which the average path 
length is short and the total number of routes large. It is, however, not 
infallible, and counterexamples can be generated—networks in which 
the process leads to neither a minimum average path length nor a 
maximum total number of routes. 





Fig. 13— The complete routing pattern to node 1. Apply backward produc- 
tion in the order: 1, 2, 3, 4, 7, 8, 6, 5. 


VII. SUMMARY AND CONCLUSIONS 


This paper discusses methods for generating loop-free directed 
routing patterns and for detecting the presence of loops in arbitrary 
patterns. The heuristic procedures suggested seem to yield useful pat- 
terns for the size network that can be considered by hand; moreover, 
they are clearly programmable, thus allowing us to deal with large 
networks. 

The procedures and theorems we present are not addressed to the 
problem of achieving optimum traffic handling abilities of commu- 
nication networks. They are, however, a preliminary step to such ex- 
aminations and, hopefully, present an orderly and useful way of 
looking at the process of routing as it is currently practiced. These 
theorems and procedures suggest ways of modifying present routing 
practices which may be fruitful to explore. 


REFERENCES 


1. Weber, J. H., unpublished work. 

2. Hakimi, S. L., “On the Degrees of the Vertices of a Directed Graph,” J. 
Franklin Inst., 279, No. 4 (April 1965), pp. 290-308. 

3. Wernander, M. A., “Systems Engineering for Communications Networks,” 
talk at IEEE summer general meeting and nuclear radiation effects con- 
ference, Toronto, Ont., Canada, June 16-21, 1963. 

4. Bene’, V. E., unpublished work. 

5. Benes, V. E., “Programming and Control Problems Arising from Optimal 
Routing in Telephone Networks,” talk at the First International Confer- 
ee on Programming and Control, USAF Academy, Colorado, April 15-16, 
1965. 


Diode and Transistor Self-Analogues 
for Circuit Analysis 


By BERNARD T. MURPHY 
(Manuscript received September 25, 1967) 


A new method of circuit analysis based on the time-scaling of actual 
circuits has been proposed. Audio-frequency self-analogues of microwave 
frequency transistors can be constructed using charge control theory, and 
these accurately model transistor performance in the active region. Scaling 
of storage times in diodes and transistors requires multiple-lump modeling. 
The multiple lump model developed by Linvill is reformulated here on 
the basis of an analogy between charge density in the semiconductor and 
charge density in the model, rather than between carrier density and voltage. 
Only two parameters, time constants corresponding to lifetime and a 
diffusion transit time in the semiconductor, need be specified in the re- 
formulated model. This simplified multiple-lump model should be generally 
useful for device calculations. We describe a diode self-analogue which 
ts an exact physical realization of the multiple-lump model. Separation 
of active and saturation region stored charges can be achieved in a transistor 
self-analogue, so that a single-lump model can be used for the active, and 
a multiple-lump model for the saturation region stored charges. 


I. INTRODUCTION 


A new approach to circuit analysis has been poposed which allows 
high-frequency circuits to be characterized using simple low-fre- 
quency models. With this approach, nanosecond diodes and transis- 
tors can be slowed down to audio frequencies and interconnected in 
audio frequency breadboard models of the high-frequency circuits. 
Thus, high-frequency circuits can be evaluated and optimized with 
the convenience afforded by low-frequency breadboard techniques. 

According to charge-control theory,? the frequency and transient 
responses of diodes and transistors are determined by the charges 
stored within the devices. Nanosecond diodes and transistors can be 
slowed down to millisecond models simply by multiplying their stored 
charge by a factor of 10° or some other convenient value. Charge 
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storage in the devices can be classified broadly as terminal voltage 
dependent (fixed charge in depletion layers) and terminal current 
dependent (charge in transit). The former can be multiplied by con- 
necting large capacitors between the device terminals, as described 
by Levine.* The latter can be multiplied by using small resistors as 
current sensors in series with device terminals, and using the voltage 
developed across these resistors, suitably amplified, to cause charge 
storage in capacitors connected to the device.t Models thus con- 
structed have given excellent results for transistors operated in their 
active regions.* 

The models also give time-scaled storage times when representing 
diodes or transistors operated in their saturation regions, but the 
values may be in error by a factor of two or more. One difficulty is 
that the charge-control model does not provide any means for rep- 
resenting the distribution of charge throughout bulk semiconductor 
regions. It is shown here that by replacing the storage capacitors in 
the model by resistance-capacitance networks, an exact physical 
realization of the multiple-lump Linvill model can be obtained. A 
second difficulty in the case of the transistor is that the time con- 
stants for charge storage in the saturation region and in the active 
region are different in general. Section 3 describes means for over- 
coming this difficulty. 


II, DIODE SELF-ANALOGUES 


2.1 Charge-Conirol Self-Analogue 


Fig. 1 shows a simple self-analogue of a diode. The diode itself, D, 
is its own de model. Capacitor Cp is used to multiply depletion layer 


Cp 





Fig. 1 — Diode self-analogue with scaling of charge storage. 
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stored charge in D, with, 


CG, = KC; (1) 


in which C, is the junction capacitance of the diode. Cp would ideally 
have the same voltage dependence as C;. (A large area junction might 
possibly be used, with a battery to avoid forward bias). 

Resistor Rg, amplifier A (whose gain is R2/R;) and capacitor Cy, 
are used to multiply minority carrier storage effects. Minority carrier 
charge storage Qj, within the diode itself is given by 


Our =I FT (2) 


in which J, is the forward current in the diode and r is an effective 
lifetime which will usually be dominated by the bulk minority carrier 
lifetimes in the P and the N regions. The charge Q4, stored in Cy; in 
the model is 


PR. 
Or =1 Ral = NC (3) 
1 
in which parameters in the model are chosen to give 
re(™)c, Sy: (4) 
Ry 


K is the desired time-scaling factor. 
During turn-on and turn-off transients, Qj, in the diode and Q4, in 
the model obey the charge control equations 


dQu _ , _ Qy 
dt f T (5a) 
dQ, Oy 
di Kr (5b) 


in which in which J is the terminal current. In the model, a current Ip 
which is proportional to Q4,, flows in the actual diode at all times, 
maintaining the correct bias voltage on the diode at all times (assuming 
that series resistance gives negligible voltage drops). 

Provided that amplifier A has good common-mode rejection, the 
diode voltage has negligible influence on the charge stored on Cy,. 
Fig. 2 shows a slightly modified version of the analogue in which 
the full diode voltage appears across the plates of Cj,. In this modi- 
fied version stored charge on Cy; is the analogue of both the depletion 
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Rs 





Fig. 2— Alternative diode self-analogue with scaling of charge storage. 


layer stored charge and the minority carrier stored charge; therefore, 
capacitor Cp has been eliminated. This version cannot be used if non- 
linear depletion layer effects are to be represented. 

The analogue of Fig. 1 is satisfactory whenever the charge control 
equation (suitably time-scaled) satisfactorily describes the transient 
behavior of the actual diode. This is true in at least two important 
types of diode: epitaxial diodes with epitaxial layers which are thin 
compared with a diffusion length, and diodes formed in integrated 
circuits using the emitter-base junction of a transistor whose base- 
collector junction is shorted. In the second type of diode, minority 
carrier storage is confined to the thin base layer. 


2.2 Multiple Lump Analogues 


Although thin epitaxial layers are generally used for high speed 
switching diodes, such diodes are often so heavily gold doped that 
the diffusion length for minority carriers is even less than, or at 
least comparable with, the epitaxial layer thickness. In this case dif- 
fusion delays comparable with diode storage times occur during turn- 
off. The charge-control equation (5) is not satisfactory then, and 
the simple model shown in Fig. 1 is inadequate. For example, in the 
extreme case of a diode formed on a uniformly doped substrate, and 
for a reverse current equal to the forward current, Kingston’s anal- 
ysis,> which includes diffusion delays, gives a storage time tf; = 0.25 
7, and a fall time t; = 0.6 +, whereas equation (5) gives t, = 0.7 7 and 
ty = 0. 

Diffusion delays can be taken into account by using a multiple- 
lump Linvill model.® Figure 3 shows a diode self-analogue which can 
be made an exact physical realization of such a model. This self- 
analogue is justified later by its node equations which are expressed 
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in terms of the stored charge at each node, rather than the voltage: 


_ a pt Bi _ Qa -— Q dQ, 
Nodel, J—I,y=I1 Gt per 


ti Qx-1 = Qx-1 — Ox _ Ox — Qrai Oxa1 ue 
Node X. RC RC spa ae a 


Qy-1 = Qy oe Qy_ 

Node N. RC = i ee CT RC 

In this case the output of the amplifier is double-ended; each out- 
put has the polarity of the input opposite to which it is drawn. 

The Linvill model is based on an analogy between carrier density 
in the diode and voltage in an r-g-c line. The well-known continuity 
and current equations for the uniformly-doped semiconductor region 
adjacent to the transition region in a diode are 


a(N — Ne) = IN = Ne —dl/gaA 


Ot T Ox (6) 
' a a(N — Ne) 
qa - Oa (7) 
in which 
N — Ne = carrier density in excess over equilibrium density 
7 = lifetime 
Z = current 
q = particle charge 
A = cross sectional area of diode 
D = diffusion constant 
Drift current is assumed negligible. 
N 
Rs 
Zs 





Fig. 8 — Physical realization of Linvill model as a diode self-analogue. 
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The voltage-current equations for the analogous r-g-c transmission 
line are 


OF ary — wv AV 
ax = G’/V ~C ry (6a) 
aV ; 
ax = —R’'I (7a) 
in which 
Analogue of 
V = voltage N — Ne 
I = current I 
G’ = conductance (combinance) per unit length gA/r 
C’ = capacitance (storance) per unit length gA 
R’ = resistance (1/diffusance) per unit length. 1/qgAD 


The analogy can be expressed in a simpler, dimensionally-con- 
sistent way if the equations are written in terms of charge per unit 
length, Q’, in both cases. Then equations 6 and 7 become 

a(Q’ — Q’e r= Qet ol 
Y=) _ _O-Q@ at en 


Ot T ox 


_ _ 7 AQ’ — Qe’) 
I= Dy (9) 


and Equations 6a and 7a become 


a Ger (Bs) 
a _ _ aig 
ae SRC (9a) 


Analogous quantities are now 


Diode r-g-c line 
Q’ ve Qe’ Q’ 
I I 
T C’/G’ 
D 1/R’'C’ 


The length of the r-g-c line is equal to the length of the semicon- 
ductor region which it represents. Redundancy caused by introducing 
area A in Equations 6 and 7 has been removed and only two param- 
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eters of the r-g-c line need be specified. The elements 1/R, G, and C 
correspond to the diffusance, combinance, and storance in the Linvill 
model. 

Fig. 4 shows the lumped version of the r-g-c line. Its node equa- 
tions are: 


@G 


Node1, 7 = 2M 4 My Oy 


RC dt 


T > Qx-1 — Qx re Qx — Qxs dQx QxG 
Node X. + di + C 


RC RC 


4 = dQv , QvG , Qy 
Nod Qy-1 — Qv _ dQv | QvG , Qv 
eon RC ae oo. Pee 

These are expressed in terms of charge stored on the capacitors con- 
nected to each node rather than node voltages. It may be assumed 


that the line and diode have been cut into equal lengths 82, so that 
C/G =r (10) 
RC = 62°R'C' = 62° /D. (11) 


The RC product in equation 11 is the analogue of a diffusion 
transit time between sections of the diode. 

Resistor R, can be used to represent a surface with recombination 
velocity vs (or a collector junction in which carriers travel with scat- 
ter limited velocity v,). Then R, is given by 


Ox 

ib 

The node equations for the lumped network are identical with the 
node equations for the diode self-analogue represented by Fig. 3, 
provided rA = 1/G. The charge distribution in the self-analogue is 
then the same as that in the Linvill model under the same boundary 





=, . (12) 


I R S R x Ron 


— { 


Fig. 4— Multiple-lump model, 
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conditions. If the self-analogue is to operate at low frequencies, ap- 
propriate scaling factors are needed. We discuss an example in the 
appendix and in Section 2.3. It remains to be shown that the correct 
diode voltage is maintained at all times. 

In the actual diode, the applied junction voltage is 


_ AP | 
Vz = ; n Oc (13) 
in which 

Q, = charge in lump closest to the junction 

Q,. = equilibrium charge in that lump. 


In the self-analogue, neglecting internal resistance, and assuming 
that r is chosen to be small, the analogue diode voltage is: 


kT I I 
Vo5= 7 [Zeck te oo i|. 


With IJ, being the current so designated in Fig. 3, 


rA 
kT : -(1 + ra) 
paces In Ce ee eee aes 


Vea at @ i +1 
rA\ Q, 
= - In Cee et | (14) 


Bearing in mind that Q in the model is the analogue of Q — Qe in 
the diode, (13) and (14) are identical if 


ACT vat we 
Re = rA | (AGL = 1): (15) 


Rather than evaluate (15) directly, it is easier to proceed as follows. 
If Ro is chosen correctly for one condition, (13) and (14) show that 
the diode voltage will be correct under all conditions. Under de con- 
ditions it is required that the full current J should flow in the diode. 
The current through Rc should therefore replace that lost through 
the de resistance of the rGC network, and Re should be set equal to 
this de resistance. 


2.3 Numerical Example 
Figure 5b shows a two-lump self-analogue of the typical diode 
shown in Figure 5a. Numerical values for the parameters are derived 
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AFTER SCHARFETTER 
—--—— SINGLE LUMP 


-———- DOUBLE LUMP 


--—— DOUBLE LUMP, 
CORRECTED 


DIODE VOLTAGE IN VOLTS 





TIME ——> 


Fig. 5— (a) Typical epitaxial diode. (b) Model of epitaxial diode. (ce) Diode 
circuit and response curves. 


in the appendix, where the behavior of the two-lump analogue under 
transient conditions also is described. Figure 5c shows the behavior 
of the actual diode, a single-lump and a double-lump model in the 
circuit shown inset. This circuit gives equal forward and initial re- 
verse currents, and 5 volts reverse bias when the diode is switched off. 

The solution for the diode was provided by D. L. Scharfetter from 
an exact solution of the semiconductor equations using the well- 
tested procedures that he and Gummel developed.’ The single-lump 
model gives only a rough approximation to the actual diode and can- 
not be adjusted to give reasonable agreement because of its incor- 
rect storage time. The two-lump model, as first calculated from the 
diode parameters is still not in good agreement, but a 20 per cent 
reduction in the assigned value of the junction capacitance gives ex- 
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cellent results. This agreement is, in fact, better than might be ex- 
pected and results from compensating errors. 

Transient: calculations for the three-lump model give a storage time 
of 0.447. For the infinite lump model the storage time is presumably 
shorter. However, current-dependent stored charge in the “depletion” 
region was found by Scharfetter to be a significant proportion of 
the total. Also, the “depletion” layer capacitance is very large in the 
storage regions. Both effects lead to a longer storage time, and com- 
pensate for the over-estimation resulting from representing the epi- 
taxial layer by only two lumps. 


III. TRANSISTOR SELF-ANALOGUES 


3.1 Transistor Operated in Active Region 

Figure 6 shows the simplest way of time-scaling a transistor. 
Charge stored on Cy; is the analogue of control charge stored in the 
transistor. Capacitors connected between B and C, and between B 
and £ can obviously be used to represent fixed depletion layer ca- 
pacitance. They have been omitted from the diagram for the sake of 
simplicity and will not be discussed further. 

The analogy holds even under base-widening conditions® and in 
saturation provided that the control charge recombines everywhere 
within the transistor according to a single lifetime. In that case, 


Iz = Q (16) 


TQ 
in which 
IT, = base current 


Q@ = in-transit control charge 
Tq = lifetime. 


I 


I 


In the analogue r, Ry, Ro and Cy, should satisfy 


(Bc, = Kr¢q . (17) 


The distribution of controlled and control charge in high-frequency, 
double-diffused transistors is quite complicated and does not lend 
itself well to separation into “base” stored charge, “collector” stored 
charge, or even into “current” controlled charge and “voltage” con- 
trolled charge. This can be seen from the results of numerical anal- 
ysis of charge distribution in such transistors, as given by Gummel.? 
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Fig. 6 — Transistor self-analogue. 


To a first approximation however, the “current dependent” control 
charge in an npn transistor consists of (2) a component of base 
charge Qs» associated with electrons in transit through the base, and 
(2) a component of base charge Q, associated with charge in transit 
through the collector transition region. Recombination current as- 
sociated with Q, 1s very small, and equation 16 is not applicable in 
the sence described above. However, since 


Qe = I pte ? (18) 
and 
O.25 = tie (19) 
co c 2 7, BieRFE 2 
in which 
Tr = effective base lifetime 
I. = collector current 


t. = transit time through collector depletion layer, equation 16 is 
still applicable provided that 7g is interpreted as 
Tg = Ta + Se (20) 
In spite of the difficulties in modeling described above, the simple 
model illustrated in Fig. 6 has been shown to give a remarkably ac- 
curate representation of transistor operation in the active region.* 
The multiple-lump Linvill model cannot be used to represent dif- 
fusion delay in the base of the transistor in the same way that it 
was used for the diode. Suppose lump 1 represents the base section 
closest to the emitter, lump N that closest to the collector. Emitter 
junction voltage should be related to the charge stored on lump 1 if 
the effects of emitter transition region storage on high frequency 
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response is to be correctly reproduced; collector current should be 
related to charge stored on lump JN, and the two requirements con- 
flict with one another. Two-pole representation of the transistor can 
be obtained if necessary, however, using two amplifiers as Fig. 7 
shows." 


3.2 Transistor Operated in Saturation Region 

When transistors operate in their saturation regions, excess control 
charge is stored in the device. Excess control charge is also stored 
in the model because of the excess base current. But equation 16 is 
not generally valid in this region because the majority of the excess 
control charge in double-diffused transistors is stored in the collector 
region in which the lifetime ordinarily differs from that of the base. 
Also, because the collector region is much thicker than the base, a 
multiple lump model is usually needed to represent the charge dis- 
tributed throughout the collector region even though a single lump 
model is satisfactory for the base. 

If the primary problem is that the collector life-time +, is not 
equal to rz, then the model shown in Fig. 8 can be used, in which 


Rs 


h'-s oO = TO (21) 


1 


Re 


1 


R’’- ‘Cn = Te. (22) 
In this model, two time constants are obtained with a single opera- 
tional amplifier. The simplicity of the model is, however, achieved 
at the cost of loss of accuracy in the de collector voltage in saturation. 





E 


Fig. 7— Two-pole transistor self-analogue. 
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Fig. 8— Transistor self-analogue with different time constants for active and 
saturation region storages. 


As in the case of diodes, charge distribution in the collector layer 
of epitaxial transistors in saturation is likely to be troublesome. For 
example, an npn transistor which gives a storage time of 5 ns (meas- 
ured with Ipr = Iprp = Ic) has an effective lifetime of 7 ns accord- 
ing to charge control theory. This, however, gives a diffusion length 
of about 3 ym (with D = 10 cm?/sec), which is about equal to typical 
epitaxial layer thicknesses, and is inconsistent with the assumption 
of charge-control theory that charge is stored close to the junctions. 
Digital programs for circuit analysis commonly use the Ebers-Moll 
model,?° which uses a similar assumption. Predicted storage times 
are too long and current fall times are too short in this situation. A 
multiple-lump model similar to that proposed for the diode is needed 
for more accurate representation of storage and fall time. In this 
case, it is necessary to represent simultaneously (2) active region 
storage with a single-lump model and (22) saturation region storage 
with a multiple-lump model. 


In the model shown in Fig. 9, active and saturation region storages 
are separated by the following means. The combination R, , A, , Ci 
represents active region storage, R, , A». and its associated R — G — C 
network represent saturation region storage. Active region storage 
which would otherwise occur because of base current flowing in R, 
is cancelled by feedback via R,. Thus, the feedback current J, in R, , 
which flows only in Rj because of the ground connection of amplifier 
A,, is 


TRA, 
Il, = (Hess ) (23) 


in which J, = emitter current. 
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Gn Cn 





Rn 
Fig. 9 — Separation of active and saturation region storage. 


Diode D, is chosen to compensate for the forward bias voltage across 
the base-emitter junction of the transistor, and voltages across R,, R: , 
and Ri are designed to be negligible. Then if R,/R,A, = hpz + 1 
the desired cancellation will be achieved. Cancellation can only be 
partial since hyz depends on 7, . 


CONCLUSIONS 


A new technique of circuit analysis, by which high-frequency 
circuits can be evaluated and optimized using simple audio frequency 
breadboard techniques has been demonstrated. It is based on the 
use of audio frequency self-analogues of diodes and transistors, which 
can be formed by multiplying the stored charge in the devices by 
some suitable factor such as 10.* Self-analogues based on charge-con- 
trol models can satisfactorily represent (1) the base region of a tran- 
sistor and (7) diodes and transistors formed in epitaxial layers which 
are thin compared with a diffusion length. The latter condition is not 
satisfied by most switching diodes and transistors. Self-analogues can 
be constructed which are exact physical realizations of the multiple- 
lump Linvill model. These can be used to represent diodes and tran- 
sistors formed on epitaxial layers of arbitrary thickness. 


APPENDIX 


Specific Design Example 


Figure 5a shows an epitaxial diode with typical dimensions. 
Other parameters used in this section are: 


Lifetime, 7 = 3ns 
Diffusion constant, D = 10 cm?/sec 
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Epitaxial layer doping = 10*7/cm* 
Surface concentration of diffused layer = 10*°/cm?. 


The effective recombination velocity at the epitaxial interface can 
be expected to be low because outdiffusion from the substrate creates 
a built-in field which keeps minority carriers away from the inter- 
face. 1000 cm per second is used as an illustration. 

Figure 5b shows the two-lump self-analogue of the diode. We will 
assume that each of the lumps represents half of the epitaxial layer, 
so that Cy = C2. With a scaling factor of 10°: 

C 


@ = Kr = 10° X3X 10° = 3ms 


6X 2 


Ran ae 
D 


10 


Both r and A can be arbitrarily chosen. A gain of 100 and resistance 
of a few ohms, say 59, are convenient values to use in practice. Since 
rA = 1/G, this gives 1/G = 5000, and C, = 3ms G = 6 pF, and Rk = 
2.25ms/C; = 275. Finally, equation 12 with a scaled value for satura- 
tion velocity gives Ry = 25,0009, which is too high to have a signi- 
ficant effect on the model and will be neglected. 

Current-dependent charge storage will also occur in the p-layer 
and the depletion layer. The p-layer is typically diffused and a 
charge stored in it will lie close to the junction. Both effects could 
therefore be represented by an additional capacitor in parallel with 
C,, with rA reduced in value to y r A, in which j is the efficiency of 
injection into the n-layer. This complication was not introduced into 
the present model. 

The behavior of the two lump model was calculated assuming a 
steady-state forward-bias current J; followed instantaneously at t = 
o by a reverse bias current J,. Solutions for Q; and Qe during storage 
time t, are then: 


= 2.29 ms: 








mos cate GR ae i) T -t/t 
Q, = (Att +5 + Ine 
GRr _ GR + 2 2) 
T OGR oo) ert Fee ( GR 7 
I,r T 


_ GRr =o 
ets t. + 1) ex ( GR + 
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Storage time ts is defined as the time at which Q, goes to zero. After 
Ts 


a= 0 
Q>2 = Q.(t,) exp e Cn tt), 

The behavior of the model inset in Figure 5c was compared with 
that of an actual diode in the following way. The response of the 
diode was obtained from an exact computer solution of the semicon- 
ductor equation using the procedures developed by Gummel and 
Scharfetter.’ This solution is given in Figure 5c. The area of the diode, 
3.383 X 10-* square centimeters, corresponds to a current density of 
300A per square centimeter. 

Using the Lawrence-Warner’® curves, the average junction capaci- 
tance, defined as total charge per total voltage, in the range from 0 
to 5.8 volts is 0.153pF. Using the equations and numerical values for 
the two-lump model given previously, this leads to the double-lump 
curve in Figure 5c. A 20 per cent reduction of C; to obtain best fit led 
to the corrected double-lump curve. The figure shows the solutions 
for a single lump model for comparison. These results are discussed 
in Section 2.2. 


T 
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Nonlinear Distortion in Feedback Systems 


By J. M. HOLTZMAN 


(Manuscript received November 21, 1967) 


We give a method for determining the distortion effect of a nonlinearity 
in a feedback loop. 


I. INTRODUCTION 


Desoer gives an interesting analysis of distortion resulting from 
a nonlinearity of the form v + ev™ (m an odd integer) in a feedback 
loop. Sandberg considers virtually the same problem for nonlineari- 
ties with upper and lower bounds on the slope.? On page 2546 of his 
work, Sandberg suggests that Desoer’s result may be sharpened. 
Our purpose is to show how a small modification of Desoer’s analysis 
might give this sharpening and extend its applicability. 

Desoer’s method is to find conditions for a particular mapping to 
be a contraction in a ball. The method presented in another work 
is particularly suited to that problem and will be applied in this 
paper. The problem of distortion in nonlinear systems is also con- 
sidered in References 4 and 5 among other papers. 


II. NOTATION AND PRELIMINARIES 


The feedback loop (with unity feedback for simplicity) is assumed 
to be described by 


y = NL(r — y) (1) 


where the input r and output y are in some Banach space. L and N 
are linear and nonlinear operators, respectively, mapping the Banach 
space into itself. We need not, at this point, specify which Banach 
space we are working in. Rather, we refer the reader to Reference 2 
for details on two Banach spaces of interest for analysis of nonlinear 
feedback loops.* In particular, Reference 2 shows how to evaluate 


*TIt must be verified that the Banach space (or an appropriate subset) is 
mapped into itself by the nonlinearity. In particular, nonlinearities such as de- 
scribed by polynomials do not map Lz into itself. 
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the norm of the linear operator when it is defined by a convolution 
operation or by a transfer function (frequency response). 


III. THE PROBLEM OF DISTORTION 
Suppose that 


N(x) = « + €P(z). (2) 


Then the loop is linear if « = 0 and it is of interest to determine how 
the loop response differs for a nonzero e. This difference is called the 
distortion. On the other hand, we might consider some fixed |«| > 0 
and determine how small the input r must be in order that the dis- 
tortion is sufficiently small. This latter question assumes that P(x) 
is of an order less than x as 1-0. 

The following manipulation is convenient for this problem. From 
equations 1 and 2 we have 


y = L(r — y) + ePIL(r — y)). (3) 
If we assume that ([+L) has a bounded inverse where I is the identity 
map,* we obtain 

y= (I+ L)*Lr + (I+ L)*eP([L(r — y)). 
Then, if z is the response of the linearized loop, 
z2= (2+ L)"Er. 
And if é represents the distortion, 
ea me 
we have 
= eI + L) "Plz — Le] 
= M(é). 

We are thus interested in finding a fixed point of the operation M (é). 


In particular, how large is €? To solve this problem, we use a con- 
venient modification of the contraction mapping fixed point theorem. 


we 
| 


IV. THE CONTRACTION MAPPING THEOREM 


Let X be a complete metric space (with metric d) containing 
the closed set 2 and let F map Q into itself. F is a contraction map- 
ping if there is an a ¢ [0, 1) such that 


d[F (x), F(x’)] S ad(a, x’) (x, x’ 2 Q). 
*For conditions for the existence of this bounded inverse, see Reference 2, 
especially p. 2538, 
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The contraction mapping theorem (Reference 6, p. 627) states that 
if F is a contraction mapping then there is a unique x* e« © such that 
x* = F(x*), that is, x* is a fixed point of the operation F. Also, x* 
is the limit of a sequence {x,} where 


Casi = F(x) 


and x, is any element of Q. 

One aspect of using the above theorem is finding the appropriate 
set O mapped into itself. Often, the contraction mapping theorem is 
used when © is the whole space, that is, F is globally Lipschitzian. 
The analysis of Reference 1 may be viewed as a method of deter- 
mined a ball about the origin such that a mapping is a contraction 
in that ball. The general problem of simultaneously trying to de- 
termine a set mapped into itself such that a mapping is contraction 
on that set is discussed in Reference 3. The following simple theorem 
from Reference 3 is useful in this direction. 


Theorem: Let B be a Banach space. F maps B into itself and x « B. 
It ts assumed that 


(i) F has a derivative at alla e B 
(it) There is a nondecreasing function g such that if x « B, then 


I ’(e) | S gl} & = x II) 


(77) There is an ae [0, 1) such that 
a 
9 l-—a 


k = || P@o) — 2o|l- 
Then there is a unique x* ¢ Q such that 


z* = F(x*) 





IIA 
R 


where 


where 





l—ea 
Remarks: See chapter XVII of Reference 6 for a general discussion 
of differentiation in Banach spaces. 
It is often a straightforward matter to find an appropriate func- 
tion g as we shall see in the distortion problem of this paper. 


o={eze8, lla —~ a |[ Ss i \ 
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V. SOLUTION OF THE DISTORTION PROBLEM 


To apply the preceding theorem to the distortion problem of Sec- 
tion III, we first find M’(é), then a nondecreasing g such that 





are) |] Sol -— 1D = 9M ell) Go =O). (4) 
And with P(0) = 0 (for simplicity) , 
[| (Eo) — £0 || = |] eZ + L)'P@) |! 
Sk. (5) 
We must finally find an @ e[0.1) such that 
0(; 2 -) =e (6) 


With (J+Z)- and L both assumed to be bounded linear operators, 
we have 
|] 2’) || = || + L)'P’@ — Lé)L | (7) 
(assuming that P has a Fréchet derivative). It should be clear that 
our method of analysis is not restricted to nonlinearities described 
by functions of the form of v + ev™ as used in Reference 1. For the 
case of the space of continuous real valued functions with the sup 
norm* and 
Poe). =a" m an integer > 1 (not necessarily odd) (8) 
we have that 
P'(@) = me". (9) 
Then, using equations 7 and 9, 


} a’) IP She LIL + DD []-m ll @ — LE)" [I 


Smlef-([]T+DH' I dlel+u Liew yz 
a g\(|| — — & |)) 
=g(l]€|[)  & = 0). (10) 


Now, using equations 5 and 6, we obtain 
m|e|-\|( +L) || 
° € . —1 e m\m—1 
AL ii((le t+ Weibel ee elle ee ayy 





* What might be considered to be a disadvantage of using this space is that 
the norms of the linear operator are expressed in terms of impulse responses 
rather than frequency responses. 
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Condition (11) could be used in several ways. For a fixed e, we could 
determine how small ||z|] has to be in order for there to be an 
a © [0, 1) satisfying (11) and thus get a bound on the distortion || é ||. 
The emphasis in Reference 1 is in determining how small e should be 
with linearized outputs satisfying || z || < 1 and the distortion || & |] < 4 
in order for the method of successive approximations to converge at 
a given rate (a = +). The discussion on page 2546 of Sandberg’s article” 
assumes the following conditions (in our notation): 


m= 3 
| Z || = 100 (12) 
[|Z +L) || = 2. 
Then, (11) becomes 
600 | « (1 4 00 |e) < i. (13) 


If |e| is less than about 1/2900 (actually a little larger, then (13) 
is satisfied. Then, the distortion é satisfies 





Wells 
SJel$ + D7 Welle |r 
S$le|. (14) 


The bound obtained using equation 25 of Desoer’s article’ is |«| S 
1/(2150-2900), a substantially smaller bound. 


VI. CONCLUSION 


Notice that since we do not require the mapping to be a contrac- 
tion in the whole space, we only get uniqueness in Q, the ball of 
radius k/(1—a). However, the result may be strengthened by also 
seeking the largest contraction constant a, satisfying condition 11 of 
the theorem. Then the fixed point is also unique in the larger ball. 
On the other hand, uniqueness information might be available from 
another source (for example, a property of a differential equation). 

We notice that the existence of derivatives in the theorem may 
actually be relaxed if there is a nondecreasing function suitably 
bounding Lipschitz constants. We also mention the possibility of using 
transformations to facilitate the application of the result, 


308 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1968 


The following may be helpful in visualizing the application of the 
contraction mapping theorem.* Assume that condition ww of the 
theorem of Section IV is satisfied with equality, that is, 


a; ES) = 


The radius of the ball Q is k/(1 — a). Letting 








aa k 
~~ l—-a 
the condition is seen to be 
k 
gr) = 1 — 


which Fig. 1 shows pictorially. 


g(r): LIPSCHITZ 
CONSTANT ON BALL 











L<-l|F(@0)-atol|->| | 


Eats OF BALL CONTAINING 
| ALL ITERATES 2 4,;=F (xi) 
| AND A FIXED POINT | 


1 (er a Ss ee Se __ OF RADIUS r 
sy 
UPPER BOUND OF | 
CONVERGENCE RATE | 
oo (SMALLEST @) 
| r= \|x-x| 


k--BALL OF GUARANTEED UNIQUENESS-—3 


Fig. 1— Contraction mapping theorem. 
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Numerical Integration of Systems of Stiff 
Nonlinear Differential Equations 


By I. W. SANDBERG and H. SHICHMAN 
(Manuscript received November 29, 1967) 


In connection with the design of transistor circuits, for example, tt ts 
frequently necessary to obtain a numerical solution of a system of non- 
linear ordinary differential equations. In some cases, these equations 
possess a property that leads to intolerable computational requirements 
relative to the use of standard predictor-corrector techniques or general 
linear multipoint formulas of open type. 

Here we describe an alternative approach which has been used to solve 
some practical problems by permitting dramatic step-size increases (for 
example, a factor of 10°). The approach is developed in a way which 
provides some detailed understanding of why it is useful. 


I. INTRODUCTION 


In connection with the design of transistor circuits, for example, it 
is often necessary to obtain a numerical solution of a system of non- 
linear differential equations 


t+jf@,t)=0, t20,  [x(0) = 2] (1) 
in which x and f(z, -) are N-vector-valued functions of ¢. The sim- 


plest numerical-integration formula which can be in principle used 
for this purpose is Euler’s formula: 


Ynt1 = Yn oT hy? 5) n = 0 (2) 
in which h, a positive number, is the step size; yo = Xo; 
Yn = —fYn,nh) for n2 0; 


and y, is of course the approximation to z(nh) for n 2 1. 

It is frequently the case that f(z, -) possess a property that leads 
to computational requirements consistent with the use of (2) that are 
intolerable, To see clearly how this situation can arise suppose that 


oll 
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the solution of (1) is desired over some finite interval [0, +], and con- 
sider the very special case in which f(z, t) = Ar with A an N x N 
matrix possessing distinct eigenvalues {a,;} all of which have positive 
real parts. Then using the fact there exists a nonsingular transforma- 
tion 7 such that 


A=TDT"', D = diag (a,, a, °°: , Gy) (3) 
we have 
Yor. = Ty — hD)T"Yyn ’ n20, [Yo = Xo] (4) 
in which ly is the identity matrix of order N. From (4) 
Yn = T(1y coal hD)*T~*x ; k = 0. (5) 
Since 
x(kh) = Te?’ Ta, , k=0 . (6) 
it is evident that the numerical solution is “acceptable” if h is so small 
that (1 — ha,)* is an “acceptable” approximation to e~*" for all 7 


and all values of k for which 0 S$ kh S 7+. On the other hand if for 
some value of 7 


|1—ha;| =1, or |1—ha;|>1 


then for at least one initial condition vector x , {|] yx {]}% (||-|{ denotes 
the usual Euclidian norm) does not approach zero as k — © or is 
unbounded, respectively [that is, (2) is numerically unstable]. Therefore 
if the sequence {y,} defined by (4) is to be a good approximation to 


the samples of the solution of (1) with f(z, t) = Az, it is certainly 
necessary that 


11 —ha;{ <1 for all7. (7) 


Moreover, in order to fully determine the character of the solution of 
the differential equation, it is reasonable to assume that 7, the length 
of the interval over which the solution is desired, is proportional (by 
some factor e such as 3 or 10) to the reciprocal of min; Re(a;) (that 
is, proportional to the largest time constant of the system). Thus in 
addition to (7) we have 


r= c{min Re (a,)]°. (8) 


A lower bound on the number of evaluations of (2) necessary to 
compute the solution is +/h where h satisfies (7). If all of the a; are 
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real, the smallest lower bound is simply 


max (a;) 
2 
2° min (a;) 8) 
It is a simple matter to give examples of, for instance, positive-element 
linear RC networks governed by a state equation of the form 4 + Ax = 0 
for which the bound (9) can be made arbitrarily large by choosing the 
value of one capacitor to be arbitrarily small. Thus, from the practical 
viewpoint, computation based on (2) can be impossible as a result of 
the presence of parasitic circuit elements that have no really signifi- 
cant effect on the circuit performance! It is not surprising therefore 
that a more complex and pressing problem of the same type arises in 
connection with the numerical solution of the nonlinear differential 
equations of transistor circuits, as a result of, for example, the para- 
sitic capacitors associated with the models of transistors. For many 
practical circuits of this type, computation time estimates, based upon 
use of (2) and a modern high-speed computer, are about 1000 hours. 
The well-known basic problem described above arises not only in 
connection of the use of (2), but (as can easily be shown) is en- 
countered also in attempts to use more general integration formulas 
of open type? 2 


Ynti = > UYn-~ +h sy D.Yn—n 3 (10) 
k=0 = 
or predictor-corrector techniques? ? such as 


Yor” = Yar + 2hyf 
ees = Yn = shyt ae Goa); 


The fundamental difficulty associated with the integration of “stiff 
equations” results from the restrictions that must be imposed on h in 
order to insure numerical stability. 

The purpose of this paper is to consider the properties of alterna- 
tive numerical methods for obtaining solutions of equations of the 
form (1). Our principal objective is to present some analytical results 
that shed some light on the properties of a class of numerical-integra- 
tion techniques that have been used to solve practical transistor circuit 
problems by permitting dramatic step-size increases (for example, a 
factor of 10*) relative to the methods defined by (10) and (11). 

More explicitly, attention is focused on “large-h algorithms” based 


(11) 
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on, or derived from, the standard formula of closed type 


Yntr = Yn + Ayres (12) 


which is a special case of the general multipoint formula 


Yat = do QYn—-~ h .3 b. Yh. (13) 
=0 k=-1 


with b_, ¥ 0. There is an extensive body of information concerning 
(12) in the numerical-analysis literature only for the case in which 
his “sufficiently small.” 


II. INTEGRATION FORMULA 


If we use the numerical-integration formula 


Yntr = Un + hynay (12) 


- in an attempt to compute the solution of (1), then y,41 is defined 
implicitly in terms of y, through 


Ynt1 + hf[ynss , (n + 1)h] = Un » n = 0, [Yo = al; (13) 

For the special case considered in Section I, in which f(x, t) = Ax 
and A = TDT-, we have 

You = Ty + hD) Ty, , n=z0 (14) 


and to the extent that (1 + ha,)~* is a good approximation to e7*", 
(13) generates an acceptable numerical solution of the differential 
equation. More explicitly (13) generates the exact solution of the 
differential equation 


é¢+Ax=0, t20, ~~ [x(0) = a] (15) 


in which A = TDT™' and e~™ = (1y + AD)". 

Let us suppose now that all of the a; are real and that ha, is very 
small relative to unity for 7 belonging to a proper subset $ of N = 
{1, 2, --- , N}, and that ha; is very large relative to unity for 7 belong 
to the complement § of $ with respect to 9t. Then for all ze $, G; , the 
ith element of D is very nearly a; , while for all 7 ¢ 8, a; < a; and 4; 
is very much larger than all of the @; for which? ¢ 8. 

In other words, roughly speaking, (13) generates a solution to a 
differential equation governing a system similar to that governed by 
z + Ax = 0; the former system has virtually the same low-frequency 
performance and less pronounced high-frequency performance. To 
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look at the situation in still another way, in using (13) we are able 
to (z) break away from an extremely restrictive requirement on h for 
numerical stability, such as (7), and (zz) trade step-size tor accuracy 
of high-frequency solution components. 

The simple heuristic argument given above suggests that the use 
of (12) can lead to a considerable increase in permissible step sizes 
for a class of nonlinear transistor circuit problems in which typically 
the Jacobian matrix df(x, t)/dx of f(x, t) along the solution of (1) 
possesses only real eigenvalues which are widely separated. This argu- 
ment is supported by a proposition, proved in Section IV, which is 
concerned with the case in which there exists a constant m > 0 such 
that (with ( , -) denoting the usual inner product) 


y, fy, nh) — 40, nh)) 2 m |] y II (16) 


for all m = O and all y. If this condition is satisfied for all h > 0, which 
for the scalar case is true if 


FY). 

oy on 
for all ¢ and all y, if |] f(0, ¢) || ~0ast— © or if || f(0, £) || is uniformly 
bounded on [0, «), then (as can easily be shown) |{ z(é) || ~0ast— 


or || x(£) |{ is uniformly bounded on [0, «), respectively. The Proposition 
asserts that if (16) is met and y,,, is defined for n = 0 by (13) ,then 


lA 


I] ym |] S A+ mh)™ || a || + > (1+ mh)" || Afl0, (re — k)A) || 


for all n 2 1, which implies that (13) is numerically stable for all h 
in the sense that for all h, || {(0, nh) || + 0asn—> © implies that y, — 0 
asn—> o and {|{ f(0, nh) ||}% bounded implies that {y,}% is bounded. 

Although the result stated above does not provide quantitative in- 
formation concerning the errors incurred in using (138), it does show 
under a reasonable assumption concerning f(xz,t) that unlike all for- 
mulas (10) of open type and unlike predictor-corrector methods such 
as (11), (18) defines for any step size a sequence {y,} which is con- 
sistent with either or both of two possible basic properties of the true 
solution. 

The discussion above does not take into account the fact that at 
each step errors are inevitably introduced in solving the equation 


Yn+1 + hflyn+s ? (n 7 1)h] = Yn (17) 
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for Yn41- Consider the result of using the iteration scheme 


We? =n — Af, @+ 0K, y= 4% 


which is the usual method described 1? in connection with the theory 
of integration formulas of closed type. For the linear case [that is, 
for f(a, t) = Ax], 


k 
Yer = Dy (—hA)'Yyn 
7=0 


T 3 (—hD)'Ty, . (18) 


Therefore, if 7, denotes the approximation to y, computed from Yo 
after k, iterations, and if 7 denotes the approximation to y, computed 
from %, after ky iterations and so forth, then 


Gx = TO. x heenyee ee ©,,0..T Yo 


in which 


6, , = diag (35 (- hid) aes ¥ (-has)'); 


7=0 


Since (assuming now that all of the a; are real) 
kp . 
dX (—ha,)' 


provided that ha; > 2andk, 2 1, if ha; > 2 for some, then || §, || > 
as k — o for some initial condition y, , independent of the sequence 
ky , ke , +++ . Therefore the usual iteration method will reintroduce 
the auraeeisal instability for insufficiently small h which it is our objec- 
tive to avoid.* 

Let us consider now a different and more general approach of solving 
(17) for y,,, . Assume that there exists a positive constant J such that 
f(y, nh) satisfies the Lipshitz condition 


Il fy, nh) — f(y2, nh) || S q || Yi — Y2 | 


for all n 2 0 and all y, and y2 . Suppose also that the smallest eigen- 
value of the symmetric part of the Jacobian matrix of(y, nh)/dy of 





(19) 


* Similar instability results for the nonlinear case can be proved. But since 
this is hardly surprising, we shall not consider the matter further. 


NUMERICAL INTEGRATION 517 


f(y, nh) is bounded from below by m, a positive constant, for all y 
and all n. 
Ideally, we would like to determine the sequence {y,}% defined by 


Ynor thflynsr, MtVDhl=y, n20. 


Suppose that we determine instead a sequence {%,}% such that 
Yo = Yo 


lg- tll S 


and 
yx +hfyts,@+ IA] = % 


for n > 0 in which e is an arbitrary positive constant independent of 
n. In other words, suppose that at each step the local error in solving 
for Yn+1 is at most «. Then, according to Theorem 1 (Section IV) 


n-1 
| Ge — Ym [| S e(L + ADL + hm) DY (L + hm) 
k=0 
for all n > 1, which of course implies the uniform bound 
II Gn — Yell S ed + AN)", nL (20) 
Our assumption concerning df(y,nh)/dy implies that the condition 


(y, f(y, nh) — FO, nh)) = m | y |? 


of the Proposition is met. Thus it follows from the Proposition and (20) 
that if the local error in solving for y,,, is held to within e at each step, 
then the algorithm is numerically stable for all h in the sense that 
for all h (2) {|| (0, nh) ||} bounded implies that {%,}%$ is bounded, and 
(iz) |{ {(0, nh |] + 0 as n — © implies that for any 6 > 0 there exists 
an m, such that |] %, |] S e(1 + hl)(hm)™ + 6 for alln = ny 

The combination of this stability result and the heuristic argument 
of Section I strongly suggests that the following approach should per- 
mit the use of considerably increased step sizes with acceptable accu- 
racy, for many of the ‘“widely-separated eigenvalue” problems de- 
scribed earlier. Referring to (17), solve for yni1 at each step using, 
say, the Newton-Raphson technique;* iterate until some norm of 


* After the work reported here had been completed, A. N. Willson, Jr. brought 
to our attention a preprint of a paper by R. Willoughby and several of his col- 
leagues at IBM, in which an approach of this type is suggested. The preprint 
does not contain the principal results of this paper, the material of Section IV. 


518 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1968 


the difference between the last two iterates is not greater than some 
small prescribed constant. 

In particular, notice that for f(z, t) = Az, this approach, using the 
Newton-Raphson iteration procedure, reduces to the use of the for- 
mula Ynu1 = (ly + hA)-y, (that is, to equation 14). 

The technique described above has provided a significant reduction 
in total computation time for several types of practical problems. It 
was used, for example, to solve the system of differential equations 
governing the circuit of Fig. 1, an oscillator designed to supply a 1 
ke signal. The 16 G Western Electric 100 Mc. silicon transistor of 
Fig. 1 was represented by a charge-control model (see Section 6.2, 
pp. 556-557 of Koehler*) using two nonlinear charge-controlled voltage 
sources, with the result that the system of equations for the circuit 
is of order 5. 

Motivated by the fact that the local-truncation error for formula (12) 
is 4h’£(€) for some é~ e[nh, (n + 1)h], the following method was used 
(for this problem as well as for others) to control the step size. Let e 
denote the largest of the magnitudes of the elements of the vector 
of second differences associated with the most recently computed point. 
If e « [42, é] (for this problem é@ was taken to be 10°*), then the point 
is accepted; if e > @, then the point is rejected and the calculation 
is repeated with h replaced with 3h. If e < 4é, then the point is accepted 
and h is replaced by 2h in the computation of the next point. Average 
step-size increases of about 10* (relative to, for example, the use of a 
forth-order predictor-corrector method) were obtained for this problem 
(see Fig. 2). 


3.0V 





0.01 F 0.01 pL F 





2.5V 


Fig. 1— One-kilocycle oscillator using a 16G “100 megacycle” transistor, 
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OUTPUT IN VOLTS 











TIME IN MILLISECONDS 


Fig. 2— Comparison of computed and experimental response of the oscillator 
shown in Fig. 1. 
III. AN EXPLICIT INTEGRATION FORMULA 


Of particular interest in connection with the approach described 
above is the numerical-integration formula 


Yaur = Yn — {ly + Af'lyn, (nm + 1)A]} “hflyn » (n + LAI, 
n2Zz QO, [Yo = 2X] (21) 


which is obtained from 


Yuet + hflY as ) (n + 1A] = Y, (22) 


by replacing Y, by yn and using as the approximation yn, to 
Yn41 the result obtained by using one step of the Newton-Raphson 
iteration scheme with y, the initial point. That is, with 


Yatr = Un — [Q’@) le-val QC) Lenn : 


In spite of its relative simplicity, it has been found that formula 
(21) is useful for solving problems of the type that we have been 
considering. For the problem of Fig. 1, it has led to an average step 
size increase of about 10°. 


(23) 
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In view of the simplicity of formula (21), and especially in view of 
the fact that yn41 is defined explicitly in terms of yn, it deserves 
special consideration. 

Theorem 2 (Section IV) asserts that for any h > 0 there’ exist 
positive constants k, and ke such that ky < 1 and 


lym I] SFr UL Yo || + hk a || #10, ( — &)A) || (24) 


for n = 1, provided that the Jacobian matrix df(y, nh)/dy satisfies 
certain conditions. For the scalar case, these conditions reduce to: 
_(1) there exist positive constants k and m such that 


m =< Um < ;, 
| | oy 
for all y and alln > 1 
a | 9 Ofy, nh) _ aflay, mh) . 
Gi) - eee 


for all y,n = land ae (0, 1]. 
Clearly, under these conditions, y, > 0 asn— oo if f(0, nh) > 0 as 
n— o and {yn} is bounded if {||f(0, nh)|]} is bounded. 

The function f(y, nh) of Fig. 3 is one for which conditions (7) and 
(11) are clearly met. If condition (22) is not met, then (24) need not 
follow. To show this, consider, for example, the function of Fig. 4 





‘~-SLOPE =1 


Fig. 3 — Definition of f(y, nh) for all n. 
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which meets condition (7), but not condition (a). We have from (21): 


I 


h=1 and y=1 imply that y, = —1 


and 


h=1 and y, = —1 imply that y. = 1 
from which it is clear that for this function Yn = (—1)" if h = 1 and 
Yo = 1, which of course implies [here f(0, nh) = 0 for all n] that (24) 
is not satisfied. Thus we see that if condition (72) is not met, then (24) 
need not follow. 






f(y,nh) 


“~-SLOPE =1 


Fig. 4 — Alternate definition of f(y, nh) for all n. 


IV. PROPOSITION AND THEOREMS* 
Proposition: If {yn} satisfies 


Ynsi tAflYnsr, (1 + Dh] = Yn, = 
and if there exists an m > 0 such that 


y, fy, nh) — fO, nh)) = my |P, n2Z0 
for all real y, then 


IV 
o 


Nye HS + mh) Hye I DE mb) II Hff0, (x — BHI | 
forn = 1. 


* Throughout this section, ||-|]| denotes the usual Euclidean Norm and (-,:) 
denotes the corresponding usual scalar product, 
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Proof: Clearly, 
(Yost » Yn — hf(0, (n a3 1)h)) 
= (Yur1 »Ynvr + Aflynsr , ( + Yh] — Af{O, (n + 1)h)) 
2 (1 + mh) || ys II’, 
and, by the Schwarz inequality, 
(Ynsrs Yn — Af[O, (n + 1)A]) S Hi Ynt1 Il] Yn \| 


+ ff ynsr [I] AflO, (a + 1A} I. 
Thus 


II Yea I] SL + mh) | yn |] + OL + mh)™ |] flO, (m + 1A] II 


from which we have 


ye HP SF mh)™ [Io Do EE may" |L AGL0, @ — FDA] || 
for n = 1, which completes the proof. 


Definition: Let \(y, nh) denote the smallest eigenvalue of the symmetric 
part of of(y, nh)/dy. 


Theorem 1: Suppose that there exists a constant m such that My, nh) 2 
m > 0 for alln = O and all y, and that there exists a constant | such that 


ll f(y, mh) — fly2, mh) || S Ulm — ye Il 
for alln = Oand all y, and y.. If {y,} satisfies 
Yner thflynir, M+ Dhl =y, on 


IV 
o 


af, with € a positive constant, {§,} satisfies 
lj. — yt || Se for nZ=O with 


Yeu + Ahfyk., m+ 1h) = Fp 
then 


|| Gi = Ys || < el + hm)™ l| Go — Yo || 
i Anan OO aS YE ha 
k=0 


forn 2 1, 
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Proof: We have for n = 0: 


Gnr + AALGnar + (Yer — Gui), (M + U)R} = Ga + Gner — Yrs) 
and 
Yat + fl bh Yh — Gad; @ FDA] 
= Yn + Aflynsr + (Yt — Guar), + DA] — Aflynsr , (@ + DA]. 
Therefore 
Gn+1 — Ynor + AflGnsr + (Ysr — Jnr), M+ VA] 
— Rflynsr + Yater — Gnor), ( + VAI 
= Gn — Yn + Gaer — Yer) + hfYasr , (n + 1A) 
— Aflynsr + (Yrs — Gner), + 1A). (25) 
With f% the symmetric part of dof(y, nh)/dy, the inner-product of 


(Gns1 — Yn+i) With the left side of (25) is 
1 
|| Gn+1 — Yn+i ||? + WA Gn — Yn+1 » [ fi{olGne1 + (yms1 aa Gn+1)] 


5 a oe 0) [Yns1 + (yt — Your], M+ DA} daGaar — tms)) , (26) 


since 


fGn+1 = (ym. = Gn+1); (n oe Dh] a flyn+1 + (Yiter — Gn+1); (n a 1)h] 
2 | * afly, (n + Ih) 
0 fe] 


da(Gn+1 _ Uaer)ie 
y veal 1+(1-a)f 


Expression (26) is bounded from below by 


(1 + hm) |] Gnsx — Yner [I 


By the Schwarz inequality, the inner-product of (G41 — Yn+1) with 
the right side of (25) is bounded from above by 


I Gasi Yer IIe [) Ga — He II 
+ || Yn+1 — Yn+t |-|| Te me l| a || Ynsr — Yn+1 [| 
| hflynss ’ (n at 1)h] = hflynss os (Yer a a), (n a 1)h] (I, 


which is bounded from above by 


[| Gara — Yost [Tel] Ge — Ym IP A [Gav — Ynsr |] Ce + Ale). 
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Thus, 
ll Goss — Your |S (Lt hm) | Ge — ya ll + CL + hm) (1 + hile, 
from which it follows that 


Il Gn — Yn || S A + hm)™ || Go — yo II 
+ (1+ hm) "7 + hil) > (1 + hm)™* 


for alln 2 1. 
Theorem 2: If {y,} satisfies 
Yas = Yn — {Iv + hf'lyn, (m + AN} “Aflyn , (2 + DA] 
for n = 0, tf 
(1) there exists a constant k < © such that 


of(y, nh) 
oy 














<k 


for alln = landally 

(it) there exists a constant m > 0 such that \(y, nh) = mforalln 2 1 
and all y 

(itt) with F ~ hf'[y, (n + Ih] and F, = hf'lay, (n + Dhl, the sym- 
metric part of {(2F — F.)F{} is* nonnegative definite for all y, 
all n, and all a ¢ [0, 1], 


then there exist positive constants k, and k. such that ky < 1 and 
ym I SRE [yo |] + aks DY hi |] (10, @% — kh] || forall n 21. 
Proof: We have 
Yns1 = Yn — {Iv + hf'lyn, (m+ Dh" {hflyn, (m+ Lh] — hfl0, (n + 1h]} 
— {ly + Af'lyn, (n + DA]} “AO, (rn + 1A]; 


hence 


“{] yn I 


+ || dw + FY [I] Afl0,@ + Dr] |] 7) 











| yner {I < | Iv — ly +7 [Fe da 





* The superscript t denotes matrix transposition. 
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with the understanding that F and F, are evaluated at y = yy, since 


hile » + Yh] = Wfl0, (+ DAL = fo Aiton, Ge + DK) devs. 


We now prove that there exists k, ¢ (0, 1) such that 


| 


for all n and all yp. 
From condition (7), with V an arbitrary N-vector, 


(QF — FV, FV) 20 


Sk, 











Peete ak Rs 
0 


QF'V, FLV) — (FLV, FLV) = 0 
which implies that 
LAV IPS QR VY pr ev Se AP 
or 
De ieee oe) ea a 


In view of conditions (z) and (72), it is evident that there exists a 
~ « (0, 1) such that 


—2FiV,V)s—-(1- || VIP -— 20 -—2¢F'V, V) -— 1-8 || Pv |p 
for alla, n, yn, and V. Therefore 
| FS — FV |) — FV |)? — 20F.Y, V) 
<~-(-24|[V]|P- 20 - )F'V,V)-—-8 || FV IP 
which is the same as 
VIP + LG! -— Fav IP + 2(F — FV, V) 
S El] VIP + QR'V, V) + é |] FV II 
or 
[| dw + F' -— FOV IP Sl] dy + FV IP. 
With U = (ly + F*)V, we have 
[| dy + FY ~ F2dy + F’)"U |P sé] U IP. (28) 
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Since (28) is satisfied for all U, 


| dw + F! — Fw + F')™ || S ky 
with k, = (&)?. However, 
ile Be ae Se | il Cas Clay a ae), 


and 











1 
Ly (ly + Fy f F, de 
0 


1 
sf \l\Gv+ Fy Gy+F-F,) Ilda sh. 
0 
Consider now || (1y+/']|. Since for any V 


| (dw + F)V |? = || VIP + 20PV, V) + | FV ||? = (1 + 2hm) |] V IP, 
it follows at once that 


I] (iw + FY || S (1 + 2hm)* 
Thus with k, = (1+2hm)-3 


I Ynea I] S Hx |] ym |] A he |] Rfl0, (wm + 1A] || 


from which we obtain the bound on ||y,|| stated in the theorem. 


V. FINAL REMARKS 


The algorithm described in this paper is a marriage of two stand- 
ard techniques, the use of a well-known closed-type numerical-inte- 
gration formula and the Newton-Raphson iteration procedure. It is 
clear that the approach is of use in connection with a certain class of 
practical problems, and, what is of at least as much importance, we 
have some detailed understanding of why the algorithm is useful. 

It is also clear that some natural generalizations and extensions 
of the approach, such as using different closed-type formulas* or 
different methods of solving systems of nonlinear equations, will lead 
to more efficient techniques. Finally, since there are several alternate 
approaches available which are also of use in certain cases (see Pope, 


i 


Az, we have Yau = TET yn, in whic = diag [(2 — oe Fueaey 
hay)(2 + hay)“|(T and the a; are defined in Section I). In view “of the relation 
between the local-truncation errors of the trapezoidal rule and formula (12), this 
suggests that for nonlinear problems the trapezoidal tule should permit larger step 
sizes for the same accuracy when the “fast components” of the solution have decayed 
to a very low level. 


* For example, for the trapezoidal rule yay: = Yn + $h(yn’ + pun ) and for Tet _ 
2- 
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for example)* much work directed toward the comparison of avail- 
able methods is needed. 
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An Upper Bound on the Zero-Crossing 
Distribution* 
By NICHOLAS A. STRAKHOV and LUDWIK KURZ} 


Let Q(T) equal the probability that a random process, x(t), does not 
cross the zero axis in a given interval of length T. A family of upper bounds 
for Q(T) ts derived with only weak restrictions imposed on x(t) and tt ts 
shown that for gaussian random processes only one member of the family 
provides useful formulae. Specific results are obtained for x(t) representing 
a number of interesting random processes. 


I. INTRODUCTION 

Let Q(T) equal the probability that a random process, x(t), does 
not cross the zero axis in a given interval of length 7. The problem 
of determining Q(T) (and related functions) has important appli- 
cations in communications theory and has been investigated by 
many authors.’-* Reference 5 gives an extensive bibliography of 
most of the related work on this subject prior to 1962. Despite all 
this effort, Q(7') is known only when x(t) is a simple nongaussian 
process (such as a process whose zero-crossings obey the Poisson 
distribution) or a stationary gaussian zero-mean process with one 
of four explicit correlation functions.® * Most of the rest of the results 
obtained are either approximate or form upper or lower bounds.° 

In this paper, we develop a whole family of upper bounds on 
Q(T). For computational purposes, however, only one member of 
the family has been found to provide useful results for most cases 
of interest. 


II, DERIVATION OF AN UPPER BOUND ON Q(T) 
Consider the transformation 
1 T 
2r = | sgn [a(t)] dt (1) 
0 


* An abbreviated version of this paper was presented at the Fifth Allerton 
Conference on System and Circuit Theory, Monticello, Ill., October 4, 1967. 
t New York University. 


529 


530 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1968 


where x(t) is a sample function of a stochastic process,* 7’ is a fixed 
observation interval, and zr is a random variable defined by the 
stochastic integral (1). The function sgn[x(t)] is defined as 


+1, x>0 
sen [z] = 0, z=0 
i eos 0; 
Since zr is a random variable, it has a cumulative distribution func- 


tion, P (zr), associated with it. From (1), two properties of P(zr) are 
immediately apparent, regardless of the statistics governing x(t): 


(i) Pér) =0 for z, < -1 
and (2) 
Pr) =1 for z, > 1 
(22) lim [PL + €-) — PQ — 8} = Qu) (3) 
and 
bed Ed OLS Ort) (4) 
where 
Qu(T) = Prob {x(t) 20 for OSt Ss T} (5) 
and 
Q(T) = Prob {2() $0 for O05t 8 T}. (6) 


Obviously, Q(T) as defined previously is related to the last two 
quantities by 


Q(T) = QUT) + Q1(7). (7) 
If x(t) is a symmetrict process, then 
QT) = Q:(T) = 3Q(7). (8) 


As a consequence of properties (i) and (iw), P(zr) can be rep- 
resented by 


Plér) = Ger) + Q1(T)uler + 1) + Qo(T)uler — 1) (9) 


* Throughout this paper, we assume that almost all sample functions of the 
stochastic pocess are continuous. Thus, (1), (5), and (6) are well defined. 

t The stochastic process x(t) will be ‘called symmetric if the probability meas- 
ures that govern it also govern the process —z(t). 
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where G(z,r) is continuous at z7 = +1 and 
wave i for « <0 
1 for x20. 


We assume throughout that the function G(zr) is not identically 
equal to zero. If it were, then it would be easy to show that Q(T) 
is known exactly; that is, Q(7) = 1. 

Next, consider the even-order moments of P(zr), denoted by the 
Stieltjes integral 


1 
dc] aPC) EH OA 2) han, (10) 
-1 


By substituting (9) into (10) one obtains 


an = [ e# dGGs) + QP) + O00) &=0,1,2-. Y 


Neglecting the first term in the right side of (11) (which is always 
positive) and taking (7) into account leads to a family of upper 
bounds for Q(7’) expressed by 


QT) Sq &=0,1,2,---. (12) 
For k = 0, (12) reduces to the obvious result 
Q(T) s 1. 


Before discussing the usefuless of the inequality (12), an expression 
for the moments will be derived. 
From its definition, (10), qo, can be expressed as 


dar = El\zr} 


where H{-} denotes the expected value of the quantity enclosed in 
braces. Substitution for z7 from (1) results in 


ih T T T 
da = at By [ff weute) + ved at, dt -- at} (13) 
0 0 0 
where 
y(t;) = sgn [a(t,)] a= 1, 2, --- , Qk. 


Interchanging the order of integration and expectation yields 


1 T T T 
i ze | i a : RG ite Ghd ea. (14) 
0 0 
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where 


Rilr, ta, +, bee) = EfyG)y(e) +++ yh) }. (15) 
We now make some remarks concerning the ordering of the family 
of inequalities expressed by (12). 
Denote the first term in the right side of (11) by sox, or 


1 
Sox =i, zr dG(zr). (16) 
-1 


We next establish that so, > Sorse (kK = 0,1, 2,...) which in turn 
establishes 

1>@>qa::: > Q(T). ; (17) 
The former inequality follows directly from 


1 
Sora = 22"* dG(zr) 


—-1 


IA 


/ 2" dG(zr) 
-1 


Sox 
with equality if, and only if, G(z7) is of the form 

Ger) = Auer + 1) + Bue, — 1). (18) 
Since G(z,) is continuous at 27 = -£1, equality is not possible and there- 
fore 

Son+2 < Sox 

which, together with (11), establishes (17). We next establish the 
readily proven fact that 


lim sx. < ¢, e>0O (19) 


kaw 


and therefore, 
lim gor = Q(T). 


We begin by choosing an a(Ko) > 0 such that 


-1+a(ko) € 
7 ar’ dG(zr) < 5? e> 0. 
—-1 

This can always be done because G(zr) is continuous at z7 = —1. 
Using the definition of s.,,(16), obvious symmetry properties, and 
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the fact that 


—-lta(ko) 


—l+ta(ko) 
i; 2 dGler) < Zt dGlen) 
-1 


Fab 


for each k > ko, it follows immediately that 


1 1-a(ko) 
lim [2 dG(r) < + lim 2% dG(er). 
ko aa E k-+0 ~—1+a(ko) 

Since the sequence of functions {zr}, k = 0, 1, 2, --- is uniformly 
convergent to zero on the interval [—1 + a(ko), 1 — a(k,)], the limit 
and the integral may be interchanged yielding (19). 

In light of (17) and (19), it appears that (12) should be evaluated 
for as large a value of k as possible. For the special case when x(é) 
is a stationary gaussian random process (assumed to be zero-mean 
without loss of generality), it does not seem to be possible to evaluate 
d2x for k > 1 as evidenced by the following discussion. 

As shown by McFadden,’ the quantity R(t, , t2, --: , t), defined 
in (15), is equal to the sum of some simple terms plus a quantity P,(r), 
which is defined as 


= m ss Ia 

Pi) = Qn)? [rl i dx, ++ / dx, exp | -4 Dy rite | 
(0) 0 47 

where r is a covariance matrix with elements 


ty = rt; — t;) = E{xti,)x(t,)}, 9 Ay een 
| r | is the determinant of r 


> r7)2.2; is the quadratic form associated with the inverse of r 


and 
a= 20); a = LDP ee Mn 


In other words, P,,(r) is the probability that the n jointly distributed 
gaussian random variables, x(¢;) (¢ = 1, --- , n) are all positive. 

As discussed by McFadden,’ and even more thoroughly by Slepian,* 
expressions for P,(r) have not been obtained in terms of elementary 
functions for n > 3. Because of this fact, it seems unlikely that an 
expression for (14) with k > 1 can be obtained for a general gaussian 
process, x(t). It should be pointed out, however, that an expression 
for q, has been derived’ for p(r) = exp (—| 7 |), but without first evaluat- 
ing P,(r). This result is not included because, for this correlation 
function, Q(7) is known exactly.” 
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TIT. APPLICATION TO GAUSSIAN RANDOM PROCESS 


Assume that x(é) in (1) is a stationary, zero-mean, guussian random 
process, normalized so that p(0) = 1 where p(t) = H{x(t)x(t+7) }. 
The relationship (14) will be evaluated for the case k = 1, that is, 
for 


] T T 
n= qf [ RG, dt dt (20) 
where 
R(t, , t) = E {sen x(t.) sgn x(te)}. (21) 


The latter expression has been evaluated by many authors (see page 
58 of Lawson and Uhlenbeck’s book,® for example) and the result is 


R(t, t) = 2 sin”! = i (22) 


Substituting (22) into (20) and making the obvious simplifications 
in integration results in 


1 
w= =f  —w sin [o(tw)] du, 
0 
or, in light of (12), 


or) s* fa — asin [ptrw) du. (28) 


This result has been obtained by Slepian,> who states it as Theo- 
rem 5.* Slepian’s proof, however, is long and complicated, as opposed 
to the simplicity of the proof given here. Furthermore, extensions to 
other cases can be obtained using the new method. 


Iv. APPLICATION TO SINE WAVE PLUS GAUSSIAN RANDOM PROCESS 
We now turn to applying (12) with k = 1 to the case where 
x(t) = w(t) + A cos (2rft + ¢) (24) 
where 


w(t) is a stationary, zero-mean, gaussian random process with 
normalized correlation function, p(7), 

is a random phase constant uniformly distributed on [0, 27], 
is the sine-wave frequency, and 

is the sine-wave amplitude. 


A mo. +6 


* Notice that Slepian’s P[7, r(rv)] equals one half of our Q(T). 
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As derived in the next subsection, the result obtained is 


QZ’) S gs? + q:” (25) 

where 
a =2 [ (0 — wut) du (26) 

with 
5: ps eee A’ 1 — sin @ cos anv 
Hr) == | exp {— 4 = Sint 
A’ sin 6 — cos 2rbu 
a4 D. cos’ 6 } ua 
(27) 
b = fT 


and J,(x) = modified Bessel function of the first kind, zero order. 

The expression given below for gq; is approximate, except for 7’ = k/f, 
k = 0,1, 2, --- where it is exact and consequently the function is most 
accurate in a neighborhood of these points. In addition, the accuracy 
of the approximation improves as A increases. For small A, where the 
approximation is least accurate, (26) dominates q{’ and so very little 
error results in the upper bound, (25) for all values of A. The expression 


for qs? is given by 


b-n 
[C@-n-nsma, — nsb<n44 
0 


2 
‘ 


— (b—n — 4 + v)S(v) dv 


n+3~b 


+ [ (b—n—»80) dr, aay ee eae | 
g? ~s x 4\- i _ (b —n — 4 — v)SQ) dv (28) 

+2] —2)SQ) ar, ee ee eee 

a (b —n — 1 + v)S() dv 

— | @-n-140)S@) dv, n+¥5b<n41 
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where 
n=0,1,2,-: 
ewe 
S(v) = (1 — 4v) — = G(v) * (29) 
w-Qrv : 24 : 
G(v) =2 / on™ cosz dz aes / ph nl seaer dx (30) 
and 
A’ 
Kv) = “9 008 Qnv. (31) 


While (28) appears to be a formidable equation, it turns out to be 
easily computed, partly because S(v) does not depend on b and hence 
needs only to be computed once for each value of A. The expression 
(26), on the other hand, turns out to be time-consuming to compute, 
particularly for large values of b. 


4.1 Derivation of Upper Bound, Given by (25) 
The expression for ge, (14), with x(t) specified by (24) is 
1 vi T 
e=af | RB, t) dh dt (32) 
0 0 


with R(t, te) given by (15). Notice that the expectation in this case 
ranges over the three random variables, w(t), w(t2) and ». For con- 
venience, define 


1 wv 
Rh 4) = 52 f rh) de (33) 
where 
r(t,, &) = & {sgn [w, + a] sgn [w. + a]} (34) 
and 
w;, = w(t;) 


a; = A cos (2zrft; + ¢) z= 1, 2. 


The latter expectation is with respect to w, and we only. Writing out 
(34) in terms of the definition of H{-} results in 


1 oe rao ; 
r(t,,b) = ll — DI a . i sgn [w, +- a,] sgn [w, + a,] 
ieee 2 wi + w2 — 2o(r)wiwe, 


where 7 = fe — t1. 
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Applying Price’s theorem,’° to this equation results in 


ar(ty 5 ty) 2 { a; + a3 — jeones) 
SS on Se ee Rp ee 36 

der) ail — py a — PG) oe 
Integrating (36) and applying the appropriate boundary condition 
yields 





EOP gta) 


aay 


p(t 
rt, , b) = rh , ts) [acr)=0 + i: 
or, 
r(t, , 42) = 4 erf (a,) erf (a) 


2 2 
ay, + a2 pasa} Aes (37) 


9 p(T) 1 
at =a qa—ay oP {- 21 — a’ 


where 


erf (x) = Taal er ay: 


As a result of the natural separation of (37) into the sum of two 
quantities, define 


Q2 = 93" +43” (38) 


where, by substituting (83) and (87) into (32), the terms in (88) 
may be defined as 


2 T pv pr 
gf? = 3 i | i, erf (a,) erf (a2) de dt, dt, (39) 
rT 0 0 7 


oo 2 ots ‘id ie [ [- i 1 bts 
2 wT Jo: do Coed (ony: 


2 2 
“exp - ee dady dt, dt,. (40) 

The detailed steps of simplifying (89) and (40) are relegated to 
Appendices A and B, respectively. In Appendix A we discuss the na- 
ture of the approximation made in arriving at (28). 

Before applying the results just obtained to specific situations, a 
power series representation for (82) will be given. The series may be 
derived from (39) and (40) by expanding the integrands of these 
functions in their respective Taylor series, evaluating the resulting 
terms and adding the expansions for (39) and (40) together. This 
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procedure results in 


qe = fa (1 — u) sin™ p(Tu) du 





a A’ f' (1 — u)[eos 2xfTu — p(Tu)] tip aca’). 


2 Jo V1 — pu) 
4.2 Numerical Results and Comparisons 


We evaluated the inequality (25) with the aid of a digital computer. 
Yor the first case considered, p(r) = e7'"', f = 2 Hz, 7 ranged between 0 
and 3.5 seconds and A, the sine-wave amplitude, was either 1 or 10. 
The results of these computations are plotted in Fig. 1. 

Let us first discuss the A = 10 case. The quantity g3’(7') exhibits 
a damped oscillatory behavior much like a plot of [sin (af7')/(fT)]’ 
while the quantity g{? (7) decays toward zero quite smoothly. For 























Fig. 1— Upper bound of Q(T) for zero crossings of sine wave plus gaussian 
noise. 
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large A, the former term dominates except at its zeros located at T = k/f, 
k = 1, 2, --- . The values of qg.(7’) are shown as dotted lines on Fig. 1. 
Since Q(7') is a monotonically decreasing function, its upper bound 
can be constructed by drawing horizontal lines between the local 
minima of q.(7) and the first intersection of this line with q(T) to 
the right of the minimum. This accounts for the step-like curve drawn 
in Tig. 1 representing the upper bound for Q(7’) with A = 10. 

The curve representing A = 1, while not exhibiting as fast a decay 
as the curve for A = 10, shows some interesting features. As con- 
trasted to the last case, the g¢{?(7') term dominates the g§’(7’) term 
and consequently much of the oscillatory behavior noted earlier has 


disappeared. 
Another interesting observation can be made when the A = 1 curve 
is compared with the A = 0 curve (gaussian noise alone) for which 


Q(T) is known to equal (2/7) sin™' (e~”) when p(r) = e7!"'. (See Ref- 
erence 5.) Notice that for0 < T < 0.25 the A = 1 curve lies above 
the A = 0 curve while for 0.25 < T < 0.75 the reverse is true. 

This result can be explained by recalling that 7 = 0.25 represents 
one-half the period of cos (4zt). For intervals shorter than this, the 
sine wave is not likely to cross zero and the effect is to cause fewer 
zero crossings than would be obtained if the sine wave were absent. 
Conversely, for time intervals longer than one-half the period (7 = 0.25 
in this case), the sine wave is sure to cross zero and therefore tend to 
inerease the number of zero crossings over the noise-alone case. 

As a result of this observation, it seems reasonable to conjecture 
that for 7 greater than one half the sine-wave period Q(T), for noise 
alone, also forms an upper bound to Q(T) for the sum of a sine wave 
plus noise. 

Additional calculations were made for comparison with Cobb’s 
previously reported approximate results.’ The quantity that Cobb 
derived is an approximate expression for the probability distribution 
function of zero-crossing intervals, denoted by P,(7’). Rice gives the 
relationship between Q(T) and P,(T) in Reference 4 as 


OT) = 1 — 2 + 2%» / : : Psy aba (41) 


where v = expected number of zero crossings of (24). 
As observed in Fig. 1 and 2 of Reference 1, 
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for the large sine-wave amplitudes where Cobb’s approximation is 
valid. Thus, 








2/T y 
an=i-r+ ff Pols) dsdy. (42) 
0 0 
Cobb shows in equation 52 of Reference 1 that 
eae ee _Gba oa 
P,(s) = (24a)? exp | ga (43) 
where 
_ _ 2G + ol! 
aA 
pi = p(2fT). 


The approximation (48) is only valid for « < 1. 
Substituting (43) into (42), we obtain 


Q(T) = — 2im| 0.5 epi (1=2F) | 


+ rfow[-HE20Y]- oo (8) 


As in Reference 1, set 





p(s) = 2B (45) 
A=83 (46) 
Qrf = 0.875 rad/sec. (47) 


Figure 2 compares the approximate solution based on Cobb’s re- 
sults (44), and our upper bound (25). For 2fT' < 1, the approximate 
solution is somewhat smaller than the upper bound. For 2fT > 1, the 
approximation becomes negative and therefore of little interest while 
the upper bound gradually approaches zero as T increases. 


V. EXTENSIONS TO OTHER CASES 


The specific applications discussed should not be considered ex- 
haustive. For example, the case where x(t) is the sum of a sine wave 
plus gaussian noise could easily be extended to x(t) being the sum 
of a square wave plus gaussian noise. Although the specific formulae 
may be more complex, the general result (equations 12 and 14) is 
still applicable for x(t) nonstationary or nongaussian. 
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Fig. 2— Comparison of upper bound and approximate solution of Q(T) for 
zero crossings of sine wave plus gaussian noise. 


In addition, the derivation of the general result can be modified 
slightly to obtain a useful upper bound on the conditional probability 
that x(t) does not cross the zero axis for an interval of length 7, 
given that x(t) = 0 at the start of the interval. Slepian has inten- 
sively investigated this latter probability. See Reference 5 for his 
discussion. 
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APPENDIX A 


Derivation of q§? 


We seek a simpler expression for the term 


T T wT 
gi? = = i i i erf (a,) erf (a2) dy dt, dt, (48) 
rl 0 0 -?r 
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which was first encountered in (89) and where 
a; = A cos (2rft; + ¢), i (be (49) 


Substitution of the definition of erf (x) results in 





1 Te T 
=a | | Pi, u) dh dt (50) 
aT ty) 0 
where 
Ph, w)=f f° f[ exp[-a@? + yl dedyde. (61) 
-T 0 


The first step is to notice the following three easily established 
properties of (51): 





P(t , t) = P(t, — t) = P(h — t) (52) 
P(r) = P(x +3) a BAAD Ee ez (53) 
(1 a i) = -p(t — +): (54) 
As a result of (52), (50) may be written as 
a? = =e | "(P — 1)P(x) dr. (55) 


It is a simple matter to demonstrate that, for any function, H(r), 
satisfying the requirements of (52) through (54), 


(i41)/F 
/ | HG) ar = 0. (56) 

We next introduce an approximation to (51) that preserves properties 
(52) through (54). It is important to preserve these properties because, 
as a result of (56), if they are satisfied, 3” = Ofor T = k/f,k = 1,2,---; 
consequently, an approximation satisfying (52) through (54) will be 
accurate in a vicinity of these values of 7’. In addition, the three prop- 
erties permit fast computation of (50). 

The approximation chosen is given by 


a1 as w/2 (a12 +497)? 
[O [exp t-¥@? + y) de ay = a) Oar ae 5D) 
0 0 0 0 
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for sgn [a,] = sgn [a,], and 


a, a2 r/2 (a1?+a22)3 
/ / exp [—4(a’ + y)]dudy=& -| [ e's drda (58) 
0 0 0 0 


for sgn [a,] = —sgn [a,]. 

Essentially the approximation results in deforming the region of 
integration, as shown in Fig. 3. From this figure, it may be noticed that 
(57) is in reality an upper bound while (58) is a lower bound. Of course, 
it is easy to conceive of functions that give an upper bound to (58) 
and thus result in an upper bound for gq”. However, this results in a 
loss of properties (52) through (54). 

Evaluating the integrals appearing in (57) and (58) and then sub- 
stituting into (51) yields 


PC, te) = op Dlg, ty » te)[1 — 8] do (59) 
where 


p(y, t , ) = sgn [a, a2] (60) 
and 


Pash) =] PGsh): 


Proving that (59) possesses the properties (52) through (54) only 
requires the use of elementary integration theory and will therefore 
be omitted. As a result of these properties, (59) may be written as 


P(r) _ z f p(9, r){1 23 pear iene” 0+ cos? nen) dé (61) 
where w = 2nf 


p(@, 7) = sgn [cos 6 cos wr + 4); (62) 


__— REGION OF INTEGRATION 
as FOR APPROXIMATION 


‘N REQUIRED REGION 


eae OF INTEGRATION 


\ 
\ 
\e——- RADIUS = fay? +ao* 
\ 
\ 
(6) ao x 





Tig. 3— Deformation of region of integration for approximation in Appendix A. 
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furthermore, P(r) need only be evaluated for the range 0 < + < 1/(4f). 
Values beyond this range are related to values within the range by 
(52) through (54). To evaluate (61), an explicit expression for (62) is 
required. After studying this latter equation one finds that for 


0<7< 1/2 
T < OT 
1 for 9 < 985 WT 
—1 for 5 eres, 
p(O, 7) = 7 (63) 
1 for 25 ROS er 
—1 for ae ee a ae 


with similar expressions for r falling in the ranges 


k+1 
2, * 


Le 
i as 





hs Os ais) 


However, only the expression given by (63) is needed to evaluate 
(61) in the required range. 
Substitution of (63) into (61) results in 


= rT (2r/2)—-ar (37/2) —wr 
P(r) = ii Flr, 6) do + ii F(r, 6) do 


w/2 /2 


—1/2 7/2 


= F(r, 6) dé — | F(r, 8) vo | (64) 


(9r/2)-wrt (n/2)-wr 
where 
F(r 6) = 1 -_ po SA sey Io0e" 6+ cos? (wrt6)] 
? . 
With the help of some fundamental trigonometric identities it is easy 
to show that 


F(r, 6) = 1 = gaia oO) cos (wrt+26) (65) 
where 
2 
K(r) = = COS wr. (66) 


2 
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Substituting (65) into (64) and performing obvious simplifications 
results in 


P(r) = 5 (2m — 4wr — 26°4°?G(r)} (67) 
where 
G(r) = : e KO) cose dy [ EY fe (68) 
0 0 


Tor the time being assume n S f7 < n + 3 where n = O, 1, 2, --- 
Substituting (67) into (55) gives 


i = pr a (T — 7r)S(r) dr 


where 


—A?/2 





SG) =1— 4jr = 2-60). (69) 


Using the result (56), the latter equation equals 


wed f Ges Seay: 


Setting t = + —n/f, 
2 T-(n/f) ( - 1) ( a) 
m mw A os ch 
1) 80 at 


2 T~(n/f) 
2 (Gaye 
0 


as a result of property (53). Now substitute ¢ = v/f to obtain 


= ea 


b—n 
(i =f —n—vS80) do for n<b<nt+} (70) 
0 
where b = fT. 

This equation is the same as the first part of the final result stated 
in (28). The equations defined in (69), (68), and (66) are the same 
as (29), and (80), and (81), respectively, except for a convenient 
scale change. The rest of the results stated in (28), for various 
ranges of 7’, are derived in a similar manner as (70) was, using rela- 
tions (52), (53), or (54), as required. Because only straightforward 
operations are used to obtain these results, they will not be derived. 
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APPENDIX B 


Derivation of q{ 


The first step in simplifying the expression for gj”, as defined in (40), 
is to interchange the order of integration of the two innermost integrals 
to yield 


@) “ane | a i SB eae ai. UD 
qe te: de de. eae &, by, ly) AH Al, Aly 


where 





T 2 2 
Ba, ty ) te) = il exp | - a ee Zetia | dy. (72) 


Using the definitions of a; and az, (34), and some obvious trigono- 
metric identities, it is easy to show that 


ay + az — 2ad,d2 
= A*{cos [w(t, + t) + 2¢][a — cos (wr)] + [1 — a cos (wr)]} 


where we have set w = 2af for convenience. 
Substitution of this relationship into (72) results in 


Bla, 4 , k) = exp [—JiG, 7)] 


f exp (-Jieya cee [OG Ea) FE 2elhde 73) 


where 
Ia, 1) = 4 [ ts ce8 un) | (74) 
Jala, 7) = 1 Be (75) 


Setting 6 = w(t, + ft) + 2¢ in (78) and using the periodic properties 
of the integrand, yields 


Ba, r) = 2 exp [—Ji(a, 7)] is exp [—J.(a, 7) cos 6] dé. 


This integral is recognized as an expression for the modified Bessel 
function of the first kind (see Reference 11, page 181, Equation 4). 
And so 


B(a, r) = 2nI,[J2(a, 7)] exp [—Ji(a, 7)] (76) 
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where I,(x) is the modified Bessel function of the first kind, zero 
order. 
Since (76) is only a function of 7, one may define 
9 p(t) 1 q 77 
H(r) = e : Vice Qa. ( ) 
Furthermore, it is easy to show that H(r) = H(—r). Consequently, 
(71) can be written as 





G2 =o / (1 — w) Hw) du. 
0 


By setting > = uT in (77) and by making the change of variable a 
= sin 6, (26) and (27), which are the desired results, follow. 
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Adaptive Redundancy Removal in 
Data Transmission 


By R. W. LUCKY 


This paper suggests an adaptive filter, similar to that used in automatic 
equalization, for use as a predictor in data compression systems. It dis- 
cusses some of the applications of this adaptive predictor in digital data 
transmission. In the event of redundant data input to the system the pre- 
dictor could be used to lower the transmitted power output required for a 
given error rate or to decrease the error rate while maintaining constant 
transmitted power. The action of these redundancy-removal and restoration 
systems is analyzed in simple cases involving Markov inputs. 


I. INTRODUCTION 


In the design, analysis, and testing of data transmission systems it 
is invariably assumed that the input digits are identically distributed, 
independent random variables. However, in many actual systems the 
input digits may arise from a physical source which imposes signifi- 
cant correlations in the data train. In these cases we know that the 
entropy of the source is less than when independent digits are pre- 
sented. Accordingly, we should be able to use the redundancy in the 
input message to provide, in some sense, more efficient transmission. 
For example, we could imagine the redundancy being used to de- 
crease bandwidth, to increase speed, to lower probability of error, or 
to lower average signal power. 

Redundancy removal in analog transmission systems was investigated 
in the early 1950’s by Oliver, Kretzmer, Harrison, and Elias’~*. Each 
of these papers relied on the theory of linear prediction as developed 
by Wiener in the early 1940’s.° Figure 1 shows the basic idea. It is 
assumed that the input samples are taken from a stationary time series 
{v,}. These samples are passed through a linear filter whose output 
£, at time ¢, forms a linear prediction of the sample x, based on all 
preceding samples. The prediction #, is subtracted from the actual 
sample x, and only the error e, is passed on for further processing and 
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Fig. 1— Predictive system. 


transmission. Since the portion {4,} “removed” from the input sequence 
is a deterministic function of the error sequence, no information has 
been lost and the original sequence can be reconstructed at the receiver 
by the feedback loop shown in the figure. 

The philosophy of predictive systems has been widely studied for 
its application in bandwidth compression of telemetry data and of 
television; for example, see Kortman, Davisson, and O’Neal.*° In these 
examples the error samples ¢; are quantized and transmitted by pem. 
Because of redundancy, that is, predictability, in the source data, 
fewer digits per sample (and consequently less bandwidth) are re- 
quired for transmitting the error samples than for transmitting the 
original samples for a given fidelity of reconstruction. 

One of the difficulties with these data compression systems is in 
determining the predictor filter. Although the theory of linear predic- 
tion for stationary time series is well known, the practical determina- 
tion of the statistical properties of the input data and the realization 
of the corresponding optimum filter are nearly impossible. Generally, 
an approximate average statistical description is used for the input 
data and a considerably simplified version of the optimum filter is 
constructed. Most existing compression schemes appear to use only 
linear or zero-order extrapolation of the previous sample to form the 
prediction of the succeeding sample. More complicated and adaptive 
prediction techniques have been confined to computer-processed data. 

In this paper we describe a simply-instrumented adaptive filter for 
use as a predictor. This filter uses a finite tapped delay line whose 
coefficients are continually adjusted to provide a least squares predic- 
tion of incoming data. The coefficient settings are based on the sta- 
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tistics of a finite section of the past data (the learning period). As 
the statistics of the data during this learning period change, the 
coefficients are changed to provide an updated version of the predictor 
filter. 

Although the most obvious applications of this adaptive predictor 
would be in the transmission of television or some other very redun- 
dant analog signal, we choose here to explore its application in digital 
data transmission. In the past, little attention seems to have been 
focused on the use of prediction in digital transmission. Presumably 
this is because the most effective use of prediction would be in the 
compression of the analog wave from which the digits are taken. 

However, there do exist situations in which the input digital signal 
is not under the control of the transmission systems designer. This 
occurs notably in the design of data communications equipment. 
Although it has been common practice to use redundancy in speech 
signals to ease transmission system requirements (the TASI system 
is a dramatic example), nothing similar has been attempted with 
digital data signals. There would seem to be no compelling reason 
why any redundancy in digital signals should not be taken advantage 
of, as long as the error statistics of the output data were not ad- 
versely affected by the procedure. After describing a digital redund- 
ancy removal and restoration system we shall discuss its possible 
benefits to the customer and to the transmission plant. 


II. SYSTEM DESCRIPTION 


Figure 2 shows a digital redundancy removal and restoration scheme. 
For simplicity we assume that the input digits a, are binary, although 
the technique obviously extends to multilevel transmission. The input 
sequence is passed through a shift-register transversal filter whose tap 
gains c; have been adjusted so that the filter output 4, , where 


N 
Gn = os CpOn—k ’ (1) 
k=1 


is a linear least squares prediction of a, . This prediction is subtracted 
from the actual sample a, and only the difference e, is passed to the 
modulator for transmission. Notice that, although a, is a binary variable 
taking on the values +1, both 4, and e, are-analog. Unless the digits 
a, are uncorrelated, the error samples e, will have smaller variance than 
the unit variance of the input data. Consequently, a linear modulator 
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Fig. 2— Digital redundancy removal and restoration. 


will put out less line power in transmitting the error samples than in 
transmitting the original data. . 

After demodulation at the receiver, the missing, predictable, com- 
ponent 4, must be added to the error sample e, before slicing, in order 
to recover a, . This component is obtained by a bootstrap arrangement 
wherein the detected symbols are passed through a transversal filter 
identical to that at the transmitter in order to form the predictions 4, . 
The receiver is similar in arrangement to the circuitry used in de restora- 
tion. 

There are two relatively simple ways in which this system could 
be used to improve transmission efficiency. As shown in Figure 2 the 
system lowers the average transmitted power without appreciably 
affecting the output data error rate. In this mode of operation any 
benefit from the data redundancy is used to lower the load require- 
ments on the transmission plant. If many data sets were equipped 
with such circuitry, the average power handled by the plant would 
be lowered in a statistical fashion. Some sets, transmitting entirely 
random data, would require their normal power complement. Others, 
transmitting redundant data, would require considerably less. Notice 
that this is exactly the type of effect which now takes place for voice 
transmission. 

As the input data becomes entirely redundant in the limit, the 
transmitted power goes to zero. In this case the input data consists 
of a periodic pattern. In spite of the zero-level line signal, the pat- 
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tern is reconstructed exactly at the receiver (in the absence of noise). 
Such an eventuality would alleviate the problems now encountered 
with the transmission of periodic data. These data patterns normally 
lead to tones, that is, line spectra, in the transmission channel which 
cause certain overloading and other system malfunctions. 

Currently the problem is being treated in wideband transmission 
by the introduction of digital scramblers.2 In practice the zero- 
level transmitted signal would not be a satisfactory solution to the 
tone problem since some signal strength would be required for syn- 
chronizing and timing maintenance. However, proper design of the 
system could ensure that some minimum signal strength was main- 
tained under all circumstances. For example, a nonlinear element in 
each predictor could be used to keep the predictions smaller than 
unity. As long as the same nonlinearity were used in both transmitter 
and receiver, the data signal would be reconstructed perfectly at the 
receiver. 

The other simple way to use redundancy removal to aid transmis- 
sion would be to keep the level of transmitted power constant while 
lowering the probability of error. In this case, compensating gain 
controls would be placed at the transmitter output and at the re- 
ceiver input. These controls would be adjusted to keep the transmitted 
power constant regardless of signal redundancy. During periods of 
redundancy most of the voltage presented to the slicer at the receiver 
would come via the feedback predictor and therefore would be noise- 
less (in the absence of errors). Since the small error signal transmitted 
would be greatly amplified to keep line power constant, the total noise 
presented to the slicer after complementary deamplification would be 
much smaller than in normal transmission. Consequently, the error 
rate would be diminished during periods of redundant data trans- 
mission. 

Complementary amplification and deamplification surrounding chan- 
nel noise introduction are automatically accomplished in transmission 
over compandored facilities. Normally for these channels we would 
expect that the error rate would be independent of transmitted power 
level. In the redundancy removal system, however, this mechanism 
is defeated by using the noiseless feedback in the detection process. 

There are further uses of redundancy removal in data transmission, 
but they appear to involve more complicated system arrangements. 
For example, the bit rate and bandwidth of the data signal could be 
lowered for redundant data. This could be accomplished by slicing 
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the prediction dG, to obtain a closest digital prediction and then sub- 
tracting d, from a, in digital form. The resulting error digits could 
then be processed by run-length encoding to achieve message com- 
pression. Of course we would then need a buffer to ensure a constant 
channel bit rate. We will not discuss this type of system further here. 

Thus far we have alluded to the possible benefits of redundancy 
removal in data transmission. There is also one major drawback— 
that of error propagation. Since the estimate d, at the receiver de- 
pends on the correct reception of all previous data, the compensation 
at the receiver is perfect only in the absence of errors. When an error 
occurs, the probability of error in succeeding bits tends to be larger 
and an error propagating effect occurs. Notice that this effect does 
not depend on the particular circuit configuration for its existence, 
but is a philosophical necessity in any redundancy removal operation. 
We analyze the effect of error propagation in a simple example in 
Section V. Normally we would not expect the error propagation to 
increase the entire error rate by more than a small algebraic factor. 


JII. THE ADAPTIVE PREDICTION FILTER 


In the theory of linear prediction developed by Wiener® and others 
it is assumed that the input samples a, are taken from a stationary 
time series with known covariance function R (n), where 


E{a,,a,] = R(m — n). (2) 


The power output, which is the mean square prediction error, is 


P = Efe?| = (a, — > cats) b (3) 


The coefficients c,; k = 1, ...N, which minimize this prediction 
error, can be obtained by the solution of the NV’ simultaneous equations 


ae eer nh = 1, 2, -** NN, (4) 


In case of an infinite filter (VN = oo) the coefficients c;, and the 
prediction error are given by a method involving factoring of the 
spectral density G(f) of the input process. Under proper conditions 
the prediction error P can be expressed in the form 


P = exp | i log G(f) as| (5) 
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(See Doob for the mathematical niceties of this result.2°) Notice that 
if the input symbols are independent, G(f) = 1, | f | < %, and 
P = 1. Since the input power is also unity no gain is achieved by 
the prediction process. If, on the other hand, G(f) is not flat the pre- 
diction error, P is less than unity and power is saved. 

While the mathematics of linear prediction for stationary time 
series serve as a guide to actual system performance, it is clear that 
the assumptions are philosophically inadmissible. Furthermore, since 
the data source is outside the designer’s control, it would be extremely 
unlikely that the covariance function would be known in advance. 
For these reasons, Balakrishnan" in 1961 developed a mathematical 
formulation for a learning or adaptive predictor wherein the form of 
the prediction operator was dependent solely on the past data and 
not on any assumptions of stationarity or of prior knowledge of data 
statistics. 

In Balakrishnan’s formulation that prediction operator is chosen 
as optimum at time ¢, which works best when applied at times 
tn1,..., tz. Since all past information is available, we could “try 
out” all possible prediction operators on the previous data and select — 
the operator for which 


L 
i, = De (Gaag dn—j)'W; (6) 
j=1 


is minimum. The weights w; could be used to assign a relative im- 
portance to each past trial of the predictor. 
For our finite linear predictor we have 


L 


N 2 
EL, = >>. |e. = Yeates | W; . (7) 


j=1 
In order to develop a physical implementation for this adaptive filter 
we use & motivation based on a steepest descent approach. The deriv- 
atives of the error #, with respect to the coefficients c,, are 


E L N 
Olt, =_- — > 20,| a, = a ie allies (8) 
7=1 k=1 


OCm 





OL, _ L 

ras 2 ke ae (9) 
Notice that these derivatives can be obtained by passing the product 
of sample ay», and the error voltage e, through a filter with impulse 


response {w,;}. Thus we are led to the adaptive filter configuration 
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shown in Figure 3. This configuration is entirely similar to that cur- 
rently being used for equalization?? and for echo suppression.*® 14 

When the input samples a, are digital, the circuitry of Figure 3 is 
quite simple. The delay line becomes a shift register and the multi- 
pliers become simple polarity switches. However, the circuit is not 
limited to digital applications, but could be used in such analog 
functions as telemetry or television compression systems. 

In any event, the response of the system, involving accuracy and 
settling time as well as stability, is controlled by selection of the 
smoothing filters W(o). Basically these filters must perform an aver- 
aging followed by an integration. If the data were stationary and 
the memory L sufficiently long, the result of averaging the product 
of the error and sample voltages for the m tap coefficient would 
give (see equation 8) 


N 


Ym(t) S Eldn-men] = R(m) — Dd) ()R(m — k). (10) 


k=1 


Then these voltages would be integrated for use as tap coefficients, 
so that the governing system equations would be 


Galt) = A) R(m — S o,()R(m — | for m=1,---,N. (11 


=1 


This system would be stable for all A, since the covariance matrix, 
whose nm entry is R(n-m), must be positive definite (see Davenport 
and Root’). All voltages y,(¢) would be asymptotically reduced to 


jl iy i's 
en 
x Ve x 


Fig. 3— Adaptive prediction filter. 
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zero and the filter coefficients would asymptotically approach those 
of the optimum (least squares) linear predictor of equation (4). 

For nonstationary data and realistic filters W(w) the analysis of 
the nonlinear, multidimensional control system is extremely compli- 
cated. Let us study the dynamics of the one-dimensional system 
formed by using a one-tap predictor as a guide to the behavior of 
the system. 

In order to put this analysis into proper perspective with regard 
to the system of Figure 2 we should observe that when the input data 
statistics change abruptly, both transmitter and receiver predictors 
undergo the same transients. If the predictors are identical, these 
transients cancel exactly at the receiver summer and no loss in noise 
margin is suffered. However, the statistics of the transmitted signal 
are affected by only the transmitter predictor. Therefore, the proper 
design of the adaptive predictor is crucial to obtaining desirable line 
power statistics, but not to the performance of the entire system. 


IV. THE ONE-TAP TRANSMITTER FOR BINARY DATA 


Figure 4 shows a one-tap transmitter with a binary input signal 
of the form 


s() = = ar(t — nT) 
a=] 
Hines : OSt<T. (12) 
0 elsewhere 
The transmitted voltage is given by 
e(t) = s(t) — e(t)s(¢t — T) (13) 
where 
c(t) = Aw(é)+[s¢ — T)e(é)]. (14) 
Because of the binary nature of the input s?(¢) = 1 and thus 
c(t) = Aw(@-«[s@s(t — T) — cd]. (15) 


Let m(t) = s(t)s(t—T); then the Laplace transform solution for 
C(s) is* 


*Some liberty has been taken with the shift-register starting state. 
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Now returning to equation (13) we multiply both sides by s(t—T) 
to obtain 


(16) 


e(t)s(t — T) = m(t) — c(t). (17) 
Combining equations (16) and (17) gives 
e(t)s(t — T) = m(t)«h(d) (18) 
where 
HO = ae: (19) 


The output signal itself can be written by again multiplying equation 
(18) by s(t-T) 

et) = st — T)[m@)«h()]. (20) 
Notice that the special properties of binary sequences have been used 
in arriving at this solution, so that equation (20) does not hold for 
multilevel or analog input. 

Figure 5(a) shows the mathematically equivalent transmitter 
given by equation (20) as well as its corresponding receiver. Since 
the second multiplier does not affect the transmitted power in any 
way, both transmitter and receiver can be simplified by its removal 
to result in the equivalent represented by Figure 5(b).* This final 


* The systems differ in their noise performance, however. 
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equivalent system is amazingly simple and appears to bear little 
resemblance to the initial system of Figure 4. It is interesting to 
observe that, while the initial system was termed “adaptive,” no 
one would seriously consider its equivalent in Figure 5(b) as being 
adaptive in any sense. 

Figure 5(b) has an intriguing interpretation. The input data is 
first subjected to the nonlinear operation of delay and multiplication. 
The output of the multiplier is 


m(t) = >> a,d,—(t — nT). (21) 


This voltage has a mean value given by R(1) in the stationary case. 
If the filter W(w) has been designed as a low pass filter, then the 
filter 1/[1 + AW(w)] in the equivalent circuit is a high pass filter. 
Thus the de component of m(t) is removed before transmission and 
reinserted via a de restorer at the receiver. In other words, a nonlinear 
operation on the input signal has converted the correlation into a 
spectral line which can then be removed by a time invariant linear 
filter. It would seem that some generalization of this concept should 
be possible, but as yet none has been found. 

The equivalent circuit can be used for design purposes in selecting 












AW(s) Nara 


>a D+ 


(a) 





Fig. 5— Equivalent binary one-tap systems. (a) Equivalent system. (b) 
Simplified equivalent. 
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W(w), or for calculating line power or transient response. Here are 
the results of a few straightforward examples. 


Example 1 
Simple RC filter, dotting pattern input applied at time zero: 





1 
we=rt; a-5 
ees ce nm even (22) 
(—I, n odd 


A deterministic sequence is to be transmitted. We find that the output 
of the equivalent circuit is 





e(t)s(t — T) = la TE 7a yu + i 4 i gracene (23) 


Thus the error voltage transmitted in the original circuit becomes 


e(t) = -| > (—1)'r(é -2 ]44 | u(t) + A £ -eeceont |. (24) 


The error voltage does not approach zero because of the lack of an 
integration in the smoothing filter. 





Example 2 
Simple RC filter, markov input: 


If the input is a first order Markov process the one-tap predictor 
becomes the optimum linear predictor. (We study this case more thor- 
oughly in the next section.) The covariance function of the input time 
series is taken to be 


R(n) = R'", (25) 


Since we now are dealing with a random input, our concern is with 
the transmitted power level rather than the exact waveform as in 
the previous example. The transmitted power is the same in Figures 
4 and 5b, so we use the simpler structure of the latter diagram for 
analysis. 

When the input Markov process is subjected to delay and multipli- 
cation, it can be shown that the resultant symbols (a,a,_,) have 
mean value & and are uncorrelated. The spectral density of the 
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multiplier output m(t) is given by 


sin” et 
Sirw) = R’ de) + 1 — R*)T ———- (26) 


iz) 

2 

This spectral density can be multiplied by | H(w) |? and integrated to 
give the transmitted power. The power becomes 


2 


eee: 
“(ae Ay 


1 1 Loe 
{a + Ay t [: a+ rll al + ANT }} a 
Ideally, of course, this power should be (1—R?), but the crude RC 
filter is unable to approximate this result unless the gain is high 
and the time constant (1/«) is large. 

Better results in both examples could be achieved by an improved 
selection of the filter characteristic W(w). We can see from the 
equivalent circuit that the best choice of W(w) makes 1/[1 + AW (we) | 
an efficient high pass filter with a transmission zero at o = 0. Of 
course this must be compromised with any requirement on the filter 
response time. . 

In this section we stress the use of the equivalent circuit as a 
method of analysis rather than as an implementable system. Clearly, 
if one were to build a one-tap binary predictor, the circuit of Figure 
5(b) would be preferred to that of the original system. However we 
believe that such a restricted system would not be of great practical 
interest. 

While the implementation of the simple equivalent circuit cannot 
be extended to wider application, it is hoped that the easy analysis 
of the simple system conveys some insight into the performance of 
multiloop systems. This would be particularly true if there were 
small interaction between taps on the multiloop system. Such a 
situation would occur if the covariance R(n) decreased rapidly with n. 


va P= A) 


V. ERROR PROPAGATION 


When noise is added in the transmission channel there is some 
probability of the received digits being incorrectly detected by the 
slicer. Even though the transmitted power might have been substan- 
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tially reduced by the redundancy removal, the probability of an 
initial error is identical to that of a full power system. Once an error 
has been made, however, the probability of making subsequent errors 
is increased because of the incorrect symbol being used in redundancy 
restoration. Thus, errors tend to bunch together in the received data. 
Besides increasing the average probability of error this error propaga- 
tion considerably complicates the problems of error control in the 
entire system. 

Error propagation in: de restoration circuits has been examined by 
Zador, Aaron, and Simon.?* +7 Jt appears to be a very complicated 
problem, in general, which is even more confused by the presence of 
the adaptive, pattern sensitive filters in the redundancy removal 
system we are considering here. Therefore, we shall attempt the 
analysis of only the simplest meaningful theoretical model. Both 
transmitter and receiver will have one-tap transversal filters as shown 
in Figure 4. The input data is taken to be a binary first order Markov 
process, with zero mean and covariance 


R(n) = R'"', 


The transition matrix for this process js: 














An+1 
+1 —I1 
L+ i 1= 8 
an 2 2 
an _— = 
ag |) ee 


2 2 


The ideal linear predictor for this time series is simply Gd, = Rayn4 
and the average transmitted power using this predictor is 1 — R?. 
Since the ideal predictor uses only a single tap filter, the assumption 
of single tap filters in the actual system is not particularly restrictive. 
If additional taps were used, their gains would be small and their 
effect on error propagation would not be significant. 

We will assume that noise samples &,, uncorrelated Gaussian 
random variables with zero mean and variance o?, are added to the 
transmitted symbols in the channel. We further assume that suf- 
ficient smoothing is done at the transmitter so that the tap gain may 
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be fixed at its optimum value, R. Thus the transmitted samples are 


e, = a, — Ray-1. (28) 


Now at the receiver we shall write the received symbols as B,a;. The 
parameter £;, = +1 indicates the absence (+1) or the presence (—1) 
of an error at time ¢,. If the tap gain at the receiver is denoted by 
the parameter c, the detected symbols can be written 


Bra, = sgn [a, — a(R — cBr-1) + &l- (29) 
Thus the error parameter £;, is 
B, = sgn [Ll — a.a,1(R — cBy-1) + ml (30) 


where 7, = & a; has the same statistical properties as &,. The proba- 
bility of error at time ¢, is the probability that 8, = —1, which is the 
probability that 7, is such that the term in brackets is negative. 

Now we must turn our attention to the behavior of the receiver tap 
gain c. If no errors are made, then this gain is identical to the trans- 
mitter gain and as k > o, c > R. However, because of the presence 
of errors, the receiver tap gain tends to be different from the trans- 
mitter tap gain. At time ¢, the output voltage of the multiplier at 
the receiver is 


V_ = ByQ,Bp-10,-1 — ¢. (31) 
The random variables v; are averaged to determine the movement of 
c. Notice that, since | 8.a8%-14,1 | = 1, the magnitude of ¢ cannot 


exceed unity except as a transient starting state. This eliminates any 
possibility of a runaway in c resulting from unusual error patterns. 

We assume that the action of the loop at the receiver is to reduce 
to zero the expectation of the multiplier output voltage at time 
infinity. Thus 


E{v..] = 0 = lim E[6,a,8,-1dx.-1] — Ce . (32) 
ko 


This type of final behavior would be exhibited by systems in which 
W(w) consisted of a long term averaging followed by an integration. 
The expectation of the term in brackets in equation (82) depends on 
C.. itself, so in general we end with a fairly complicated equation requiring 
a trial and error solution for c,,. By taking the limit as k — © of the 
expectation we eliminate the dependence on time and on the initial 
probability distributions for the random variables involved. 

Define a vector random variable a = (a,, 8,) taking on the four 
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possible states (+1, +1), (41, —1), (—1, +1) and (—1, —1), denoted 
by states 1 through 4, respectively. Because a, is Markov and since the 
expression for 8, in equation (30) involves only a,, d:-1, B.-1, and 
n,, we conclude that & is also Markov. The four-by-four transition 
matrix + for & has entries p;; which may be calculated from the original 
transition matrix for the input symbols a, and from equation (30) for 
the probabilities of error in various states. Table I lists these transition 
probabilities. If the 4-entry row vector w” gives the probabilities of 
&, assuming each of the four possible states, then 


wo” = oF a. (33) 
In terms of the initial state distribution 0°” 
w™ — Dp x”. (84) 


For |R| < 1 it is clear from standard Markov chain theory (see, 
for example, Reference 18) that steady-state probabilities exist for 


TABLE I— TRANSITION PROBABILITIES FOR & = (a; , Bx) 


ie) ew? 


Q(x) = Tm 


Du = Pas = (4)[1 - (t+) | 
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the transition matrix 7, that is, ®” approaches a constant vector W as 
nm —> © independent of «°°. The steady-state probabilities of the four 
possible states can be obtained by the solution of the equations given by 


Or = w. (35) 


Some algebraic manipulation yields the probabilities 


w, = Pla. = +1,8, = +1) = ae (36) 
We = P(ag = +1, 8. = —l) =} —-—w, (37) 
w; = Pag = —1,6.2 = +) = w (38) 
w, = Pa. = -1,82= -D=F-w, (39) 


where the transition probabilities p12, pis, Poo, and Pog are given in 
Table I as functions of c, R, and o. 

The expected value of the multiplier output at time infinity can 
now be written in terms of the steady-state probabilities w; and the 
transition probabilities p,;. 


Elv..] = w, [Pir — Pie — Piz + Pra) + We[Poo + P23 — Pa — Do] 


+ ws3[Ps2 + P33 — Par — Das] + Walpar + Pas — Paz — Pas] —c. (40) 
Again some algebraic manipulation yields the result 


Ril -D14— P24 — P22 — Dio] +2 [P14 — P12] +4 [pooM12—PosPr4] 
Elv,]=— 7 * a eee el, 4] 
i) 1—DPo2+Di2—PosrtPi4 cy 


The value of the tap gain at time infinity can be found by trial and 
error. A value of c is assumed, the transition probabilities are computed 
and Elv.,| is found. The value of c for which E[v,,] = 0 isc, . Notice that 
under suitable assumptions E[v,,| gives the rate of change of the coef- 
ficient c in the dynamic action of the system. 

The probability of error after the system has settled is simply the 
probability that @,, is in a state where 8,. = —1, which is simply (we + 
W4). 


Pre + Pia 
Je Ee 2 a 2 49 
: 1 — Poo + Piz — Pos + Pia ( ) 


The transition probabilities here must be computed using cy. 
Expressions (41) and (42) have been written in terms of only 

those transition probabilities which involve errors. Thus, as « > 0, 

each of the transition probabilities in (41) and (42) approaches zero, 
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c, > R, and P, — 0. Each of these probabilities can be visualized 
as the probability that the noise (zero mean, variance o”) is greater 
than the one of these four thresholds: 


(multiply by : + ) 








(multiply by : 5 ) 





0 1—-R +e 1+Rk 


Thus po4 is the smallest transition probability, while po» is the largest. 

If the transition probabilities are small, it can be seen from equa- 
tion (42) that P, is principally determined by (pie + p14), which is 
minimized by c = R. Also we notice from equation (42) that the tap 
gain c approaches R very closely for small transition probabilities. 
In general, however, c = FR will not be the best setting to minimize 
the error probability in equation (42), nor is it the setting to which 
the loop settles. Unfortunately it appears that these are not compensat- 
ing offsets. For example, in Figure 6 we have plotted P, and E[v,.] 
against c, for a case in which R = 0.4 and o = 0.4. Although neither 








Fig. 6 — Probability of error and E[vco] ys receiver tap gain ¢, 
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effect is very significant, it can be seen that the system settles (Z[v.]| 
= 0) for a value of c somewhat smaller than R, while the minimum 
error probability is obtained at a value of c somewhat larger than J. 

In all but the most severe noise conditions the approximation of 
C. = R would be satisfactory and we would have 











OT = EB lS) - (Belt) + ell) 


(43) 


But Q(1/c) is the probability of error in the original system (no 
redundancy removal). If this probability, called P.o, is small, then 
Q(1 + 2R/c) is much smaller and we have the very good approxi- 
mation 


P 
P, er ea an % 
[Pes sma iff 2k 
ey 


The factor in the denominator gives the amplification of the original 
error rate due to error propagation. Finally if R > 1/2, then Q(1 — 
2R/c) approaches unity and we get the severe dependence upon FR 


Ay 2b 

gee = R- 

The most significant aspect of the error propagation behavior of 
the circuit is that the redundancy removal and restoration system 
has impressed the statistics of the input data (Markov here) upon 
the error statistics of the output. It is clear that this philosophy 
would hold in general. In the case of highly correlated input we would 
end with highly correlated errors. The problems of error control could 
be made quite severe in this manner. 


(44) 


(45) 
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VI. EXPERIMENTAL RESULTS 


A three-tap, adaptive transmitter and a similar receiver were de- 
signed and constructed by V. G. Koll. The system was designed for 
binary data transmission so that the multipliers in Figure 3 became 
polarity switches, while the delay line took the form of a shift register. 
The filters W(s) consisted of simple RC low pass sections followed 
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by integrators, that is, 


a 


s(s + a) 


With this choice of smoothing, the steady-state error for a periodic 
input (period 3 or less here) was zero. It was in fact observed that 
during the transmission of periodic data the transmitter could be 
disconnected with no effect on the received data pattern. 

The input data for the system was obtained by passing white Gaus- 
sian noise through a variable cutoff, low pass filter. If we assume an 
ideal low pass filter, with cutoff frequency W Hz, then the autocor- 
relation function of the filter output is 





Wo) = (46) 


sin 2a 


R,(7) = oN] onWr 


(47) 
This voltage is then sampled at rate (1/7) and subjected to infinite 
clipping so as to produce the correlated input bits. Van Vleck and 
Middleton?® show that the resulting autocorrelation is 


29s =) |sin arn |. 
He em | nw 


For a filter cutoff of 1/27 Hz the data is uncorrelated. By decreasing 
the filter cutoff frequency the redundancy in the data can be increased. 

The action of the adaptive redundancy remover is shown in Figure 
7 for two different values of filter cutoff. Notice that as the redundancy 
is increased the transmitted waveform has longer periods of near zero 
voltage where predictability is good and occasional peaks where the 
predictor is “surprised.” Except for a few minor discontinuities the 
reconstructed signal before slicing at the receiver is the same as the 
original input waveform at the transmitter. The relative power saving 
as a function of filter cutoff is shown in Figure 8. 

In order to predict system performance in Gaussian noise we make 
the crude approximation that the input process is Markov with R (1) 
as given in equation (48). According to this approximation the trans- 
mitted power should be 1 — R(1)?. This value is also shown in Figure 
8 in comparison with the actual measured power output. Since the 
exact correlation function is known, the theoretical signal power output 
could be computed precisely through equation (4). However, we have 
no corresponding means of computing the degree of error propagation 


(48) 


REDUNDANCY REMOVAL 569 





Fig. 7— Transmitted and reconstructed signals. (a) Filter cutoff w7 = 04 
{little redundancy, R(1) = 0.15]. (b) Filter cutoff #7 = 0.1 [moderate redun- 
dancy, R(1) = 0.77]. 
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for the non-Markov source. The approximate curve of signal power 
in Figure 8 is shown only as a way of evaluating the Markov ap- 
proximation for later use in predicting error propagation values. 

Bandlimited white Gaussian noise was added to the transmitted 
signal, and error rates were experimentally determined by V. G. Koll 
at a number of filter cutoff (redundancy) positions. The results of 
these tests are shown in Figure 9 in curves of probability of error 
versus signal-to-noise ratio. Beside these measured curves have been 
plotted theoretically computed curves which are based on the Markov 
approximation and on the use of equation (48) for Pe. 

Although all necessary information for performance determination 
is contained in Figure 9, it is instructive to plot two additional curves 
of probability of error versus filter cutoff. These curves are shown 
in Figure 10. In one curve the transmitter and receiver gains are held 
constant so that the line power decreases according to the curve of 
Figure 8 while the probability of error increases with increasing re- 
dundancy because of the effects of error propagation. In the other 
curve of Figure 10 the transmitter and receiver gains have been 
adjusted with increasing redundancy so as to hold line power constant. 
In this case the probability of error decreases with increasing redun- 
dancy. 


SIGNAL POWER IN DECIBELS 
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Fig. 8— Signal power saving by redundancy removal. 
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Fig. 9— Performance of redundancy removal system at various values of 
normalized filter cutoff wT. 


VII. CONCLUSION 


We have advanced two main points. First we suggest the possibility 
of using an easily-implemented adaptive predictor for data compres- 
sion systems. Second, we investigated the use of this adaptive predictor 
in digital transmission. 

We have seen that the predictor can be used to increase transmission 
efficiency for redundant data either by decreasing signal power for a 
given error rate or by decreasing probability of error for a given signal 
power. Although the required circuitry for the digital application is 
quite simple, it is nearly impossible to make an economic evaluation 
of the system because of the complete lack of knowledge of the prev- 
alence and degree of redundancy in customer input data. 
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Fig. 10 — Probability of error vs filter cutoff for constant and for free S/N. 
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Group Codes for the Gaussian Channel 


By DAVID SLEPIAN 
(Manuscript received April 27, 1967) 


A class of codes for use on the Gaussian channel, called group codes, 
is defined and investigated. Roughly speaking, all words in a group code 
are on an equal footing: each has the same error probability and the same 
disposition of neighbors. A decomposition theorem shows every group code 
to be equivalent to a direct sum of certain basic group codes generated 
by real-irreducible representations of a finite group associated with the 
code. Some theorems on distances between words in group codes are demon- 
strated. The difficult problem of finding group codes with large nearest 
neighbor distance is discussed in detail. 


I. INTRODUCTION 


In a communication model first introduced by Kotel’nikov’ in 1947, 
and independently by Shannon’ in 1948, and since studied by many 
authors,*-** messages for transmission are represented by vectors in 
a Euclidean space, §,, of n dimensions called signal space. In this 
model, known as the Gaussian channel, when X is transmitted, the 
received signal is represented by a vector Z = X + Y which consists 
of the sum of the sent vector and a noise vector Y whose components 
are independent Gaussian variates with mean zero and variance o”. 
Some physical circumstances that lead to this model, as well as further 
details, can be found in Refs. 3, 10, and 18. 

An equal-energy block code of size / for use on this Gaussian channel 
is a collection of M distinct vectors X, , X., --: , Xa in signal space 
all of the same length. We shall always suppose IM = n and that the 
vectors span S, . The length of the vectors serves to define an important 
parameter S called the average power of the code through the equation 


nS = |X,)’. (1) 


The vectors of the code are called code words or code points. Their 
termini lie on the sphere of radius ~WnS centered at the origin of §, . 
Associated with each code point X; of an equal-energy block code 
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is a region ®; of signal space called a maximum likelihood region and 


defined by 
a ={x||k- Kl s/x-x,Liei} 2) 


That is, ®; is the set of all points in §, at least as close to X; as to any 
other code word. These regions are convex flat-sided cones with apex 
at the origin. The interiors of ®, and @; are disjoint for 7 # j: the 
union of the @; is all of §,,. 

The capabilities of equal energy block codes for communicating over 
the Gaussian channel are well known. If the words of a code are presented 
equally likely and independently for transmission over the channel, 
the communication rate is 


Qa 
R= , log Mf (3) 


natural units per second where a (measured in numbers per second) 
is the rate at which vector components are transmitted. The receiver 
which minimizes the average error probability®’’ operates by asserting 
that code word X; was transmitted when the received vector Z lies 
in ®; ,7 = 1, 2, --: , M. (The received vector lies in the boundary 
of some ®, with probability 0.) When X,; was transmitted the error 
probability of this best receiver is 


1 1 | 


where &/ is the complement of ®,; . The average error probability is 


1 M 
Po = ag 2 Pa. (5) 


Upper and lower bounds are known*'®'"""'"" for P, min(/, n, S), the 
smallest attainable value of P, for an equal-energy block code with 
the indicated parameters. In the limit as n — o, these bounds lead 
to the famous capacity formula C = a/2 log (1 + S/o”) whose in- 
terpretation we suppose known. For fixed finite values of M7 and n, 
however, little is known in the general case about codes for which 
P, attains its minimal value (optimal codes). The cases M = n + 1, 


n + 2, +++ , 2n have been studied in some detail.*°""* For n = 2, 
Weber™ showed that the regular M-gon is globally optimal: for M = 
n + 1,n = 2,3, --- , it has been shown” that the regular simplex 


is optimal. No other optimal codes with n > 3 are known. 
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Recently Wyner’’ has investigated the capabilities of equal-energy 
block codes when a suboptimal receiver, known as a bounded distance 
decoder, is used. Here the regions ®; of the maximum likelihood receiver 
are replaced by spheres of radius d/2 centered on the termini of the 
code vectors X;, where d is the minimum distance between any two 
words of the code. If the received vector is not in one of these spheres, 
a decoding error is assumed. Wyner established upper and lower bounds 
on the smallest error probability attainable with an equal-energy block 
code using bounded distance decoding. In the limit as n — he ob- 
tained coding theorems and a capacity analogous to the usual ones. 
For finite JJ and n, the error probability using bounded distance 
decoding is a monotone decreasing function of the minimum distance d 
between code words of an equal-energy block code. In the general case 
little is known about equal-energy block codes with largest nearest 
neighbor distance. 

Tor equal energy block codes of 1/7 vectors spanning §, two optimiza- 
tion problems thus present themselves: to find a code for which P,, 
as given by (4) and (5), is a minimum; and to find a code with largest 
nearest neighbor distance between its code words. We have made little 
progress in solving these problems. 

In this paper we investigate instead a class of equal-energy block 
codes called group codes. It is conjectured that this class includes 
solutions to the problems just mentioned for many values of M/Z and n. 
Quite apart from these questions of optimality, however, group codes 
possess an important symmetry property that makes their study of 
interest in its own right. Roughly speaking, all code words in a group 
code are on an equal footing. This notion is made precise in the next 
section. 

Most codes that have been investigated for the Gaussian channel 
are group codes: it is likely that any code used in practice will be of 
this type. Group codes for the Gaussian channel are a natural exten- 
sion of the group codes introduced for the binary channel in Ref. 21, 
and these latter codes are obtained as a special case of the codes de- 
scribed here. 

In what follows, we define equivalence for group codes, then in- 
vestigate the possible classes of group codes. Here the theory of group 
representations plays a key role.** The appendix gives a summary 
of results needed from this field. The problem of constructing group 
codes is considered and an optimization problem of some difficulty 
is encountered. A number of interesting properties of group codes are 
disclosed. . 
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Many of the results reported here are contained in the author’s 
Bell Telephone Laboratories report of May 7, 1951, a document that 
received a limited circulation outside the Laboratories. A number of 
these results were recently rediscovered independently by J. G. Dunn 
and appear in his thesis” along with extensions in directions different 
from those reported here. The discovery of an easy decoding algorithm 
for certain group codes” has led to a revival of the author’s interest 
in this subject, and so the present paper, while in part very old, is a 
report on research now in progress. It examines the general structure 
of group codes. In a later paper we hope to give a detailed treatment 
of some group codes associated with the symmetric group. 


II. GROUP CODES 


In studying the geometric properties of equal-energy block codes, 
it is convenient to deal only with code vectors of unit length. That is, 
we set S in equation 1 equal to 1/n, and deal with normalized codes. 
To compute error probabilities associated with the use of the code, 
one must scale up the vectors by a factor ~/n8S. 

Let X, , X2, --- , Xa, be the (unit) vectors of an equal-energy block 
code. It is clear from the definition of the regions ®; and from (4) 
and (5) that P, is invariant under a rotation of the code as a whole. 
That is, if O is an arbitrary n X n orthogonal matrix and 


Xi = OX,, a= 1,2,---,M, (6) 


the error probability P/ for the code X/, --- , X4, is the same as that 
for the code X, , X2, --+ Xs, . The set of interword distances for the two 
codes is the same, and in particular both codes have the same minimum 
nearest neighbor distance d. Two codes whose vectors (with possible 
renumbering) can be related as in equation (6) are called equivalent. 
Equivalent codes have the same communication capabilities. 

We now examine in what sense the words of an equal-energy block 
code in §&, might be “alike”. Given the J7 unit vectors X; that define 
the code, the real orthogonal n by n matrix O is said to leave the code 
invariant if the Y,; are a permutation of the X; where Y; = OX;,7 = 
1, 2, --- , WM. The collection 6 = {0,, 02, --- ,0O,} of all real orthogonal 
n by n matrices that leave the code invariant clearly forms a finite* 
group under ordinary matrix multiplication. Now transformation by 


* By hypothesis, the X; span S,. An n X n orthogonal matrix is completely 
determined by its effect on a set of n vectors that span its carrier space. Since the 
words of the code are permuted along themselves by each element of 6, g S$ M!. 


GROUP CODES 579 


an orthogonal matrix preserves distances between points, so that a 
possible definition of ‘‘alikeness” for the points of the code is to require 
that in the group 6 there be elements O, , O2, --- , Oa, that transform 
any particular word, say X;, into each of the 17 vectors of the code. 
A collection of Jf unit vectors spanning &, that satisfies this condition 
will be called a growp code and denoted by the symbol {M/, n}. In a 
group code, if O; sends X, into X; and O; sends X, into X, , then 0,07" 
sends X, into X; . We have then 


Proposition 1: For a group code, the set of distances from X;, to all 
other points of the code is the same as the set of distances from X; to all 
other points of the code, 7,7 = 1, 2,--- , M. 


Each point has the same number of nearest neighbors, the same number 
of next nearest neighbors, and so on. 

The maximum likelihood regions ®; for a code are defined by equa- 
tion (2) in terms of distances from code points. Since orthogonal matrices 
leave distances invariant, it follows that for a group code a matrix 
O « 6 that sends X; into X; also sends ®, into ®;. From this fact and 
the form of (4) we have 

Epona 2: For a group cede {M, n} he maximum likelthood regions 
Ri, Ro, *** , Ra are all congruent and all words have the same error 
probabilzty, that is, Pa = Pe =e = Poy = P.. 


IfI. GENERATION AND CLASSIFICATION OF GROUP CODES 


To each matrix O of the group @ of orthogonal matrices that leaves 
a group code {J/, n} invariant, there corresponds a permutation on 
M letters, namely the permutation effected by O on the M vectors 
of the code. That these permutations form a transitive permutation 
group follows from the definition of a group code. No two different 
elements of 6 can effect the same permutation of the words of {J/, n} 
since the effect of an n X n matrix on a set of vectors spanning 8, 
completely determines the matrix. We have then 


Proposition 8: The group @ of all orthogonal n X n matrices leaving 
a group code {|M, n} invariant forms a faithful representation of (1s simply 
tsomorphic to) a transitive permutation group on M letters. 


Group codes {1/, n} do not exist for every Jf and n. For example, 
it is not hard to prove that it is impossible to arrange 5 points on the 
sphere in 3 dimensions to form a group code. Necessary and sufficient 
conditions on M and n for the existence of an {J/7, n} are not known, 
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Group codes do exist in great abundance, however, and we shall 
give examples later. Indeed, from any set of n X n orthogonal matrices 
O,,02, +--+ , Ox, that form a finite group G under matrix multiplication 
we can form a group code by choosing a unit n-vector X and forming 
the set of vectors 


X,=0.%, ¢=1,2,---,M. (7) 


Elements of G leave this configuration of vectors invariant by the 
group property. Since SG must contain the n X n unit matrix, X is 
among the collection of vectors and it is sent into each of the other 
vectors. A group code therefore results. This code may not have Jf 
distinct vectors, however, and it may not span S, . The code depends 
on the initial vector X. 

If the code has fewer than M vectors, then for some 7 # j, X; = X; 
or O;X = O;X, or O7°0;X = O,X = X for some O, « G. That is, X must 
be an eigenvector with eigenvalue unity for at least one O « G different 
from the unit matrix. The set of all such O ¢ G forms the subgroup 3¢ 
of order h of G that sends X into itself. It is easy to show that by (7) 
G generates vy = M/h distinct vectors. Since the matrices of G have 
only a finite number of eigenvectors, however, it is always possible 
to choose an X so that the M vectors (7) are distinct. 

It may not be possible, however, to choose X so that the vectors 
span &,. To discuss this matter further we must recall the notion of 
real-reducibility. A finite group of (real) orthogonal matrices § = 








O,, O2, -*+ , On, is said to be real-reducible if there exists ann X n 
real orthogonal matrix O such that for? = 1, 2,--- , M 
agp =f Ae | > 
CO ( CB; 8) 


where A; is an / by / matrix, B; is ann — lby n — | matrix,0O <I <n 
and C and D are matrices all of whose elements are zero. It is assumed 
that 1 does not depend on z. A group of real orthogonal matrices that 
is not real-reducible is said to be real-irreducible. In words, a real- 
reducible collection of matrices can be simultaneously transformed to 
block diagonal form by a real orthogonal matrix: a real-irreducible 
collection cannot be so reduced.* The reduced matrix in block form 
in equation (8) is said to be the direct sum of the two square matrices 
A, and B,. 

* In the theory of group representations (see the appendix) reducibility is usually 
defined over the field of complex numbers. The definition is as above with O replaced 


by a unitary matrix. We shall speak simply of ‘reducibility” in this case as opposed 
to “real-reducibility”’. 
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It is easy to show that if the matrices O; of equation (7) are real- 
irreducible, then the code they generate spans §, for all choices of X: 
if they are real-reducible, for some choices of X the code will not span §, . 

These comments lead to 


Proposition 4: Every real-irreducible group G = O,, O2, ++: , On 
of real orthogonal n X n matrices serves by means of equation (7) to generate 
a group code {M’, n} for each unit vector K in &, . Here M’ S M. If 
M’ < M, tt ts a divisor of M. 


Propositions 3 and 4, together with the theory of group representa- 
tions} suggest a means of classifying and generating all group codes. 
From Proposition 3 we can associate with a given group code {JZ, n} 
a unique abstract group and a faithful representation 6 of this group 
by orthogonal matrices. The code can be thought of as generated from 
one of its vectors, X, say, by the operation of the matrices of this 
representation in the manner of equation (7). Now the representation 
6 will in general be real-reducible. There will exist then a real orthogonal 
matrix O that will exhibit 6 in block form (8) as the direct sum of a 
number of real-irreducible representations. Denote this new reduced 
representation by 6’. It is easily seen that the matrices of 6’ operating 
on the vector Y = OX generate a group code Y, , Y2, °°: , Ya, equivalent 
to the originally given {17, n}. We can further regard Y as the sum 
of its projections Y*, Y’, --- on the various invariant subspaces of 
6’ indicated by its block structure. 

By the procedure just outlined, for each equivalence class of group 
codes we arrive at a particular set of real-irreducible representations, 
say 0, 02, °°: , 0; of an abstract group, each with a corresponding 
associated vector Y’, Y°, --- , Y’. We regard Y* as lying in the carrier 
space of 0; , so that if 6; is of dimension 1; , then Y* is a vector of J; 
components, i = 1, 2, --- , j. Let the length of Y* be \;. We have 
>> = 1. The 6; are determined by {M, n} only up to equivalence 
in the sense of representation theory, owing to the possibility of reduc- 
tion of 6 by different matrices O. The vectors Y* inherit some additional 
freedom owing to the M possible choices of X in the preceding paragraph. 

We can think of the {1/, n} as decomposed by the above process into 
an equivalent direct sum of 7 group codes, the zth code being generated 
by the matrices of 6; operating on the initial unit vector Z’ = Y*/), , 
t = 1, 2, --- , 7. The constituent group codes are weighted by the 
numbers A; , 42, °** , Ay in forming the direct sum code {M, n}. Notice 


+ Knowledge of the material in the appendix is necessary for understanding 
much of the remainder of this paper. 
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that some of the constituent codes may have fewer than M distinct 
words. 

Conversely, within equivalence we can construct any group code 
as the weighted direct sum of codes generated by real-irreducible groups 
of matrices after the manner of Proposition 4. In synthesizing codes 
in this manner, we may, of course, arrive at equivalent codes by several 
different constructions. The group § of Proposition 4 may be only a 
subgroup of the group of all orthogonal transformations that leave 
the code generated by G invariant. Different initial vectors operated 
on by the same group of matrices may give rise to equivalent codes. 

Every group possesses the trivial real-irreducible one-dimensional 
identity representation in which each group element is represented by 
the one-dimensional unit matrix. The inclusion of this identity rep- 
resentation in the constituent codes making up a direct sum code 
represents a waste of one dimension, since the code is then equivalent 
to one in which each code vector has the same first component. This 
first component then carries no information. By omitting the first 
component of each vector (and rescaling the length of the resultant 
vectors), a new code of dimension n — 1 is obtained with error prob- 
ability no greater than the original {M, n}. In general in what follows 
we will not be concerned with codes that contain this identity rep- 
resentation. 

We turn our attention now to the basic problem of constructing good 
group codes as the weighted sum of properly chosen group codes gen- 
erated by real-irreducible groups of orthogonal matrices. 


IV. THE INITIAL VECTOR PROBLEM AND THE FUNDAMENTAL REGION 


As in Proposition 4, let a code be constructed from a given group 
G =0,,0.,--- , Ox of orthogonal n X n matrices by means of equa- 
tion (7). We think of these matrices as a faithful representation of an 
abstract group isomorphic to the matrix group. The code obtained in 
this manner depends upon the initial X on which the matrices operate. 
The regions ®; of equation (2) and hence also P, = P,; by (4) also 
depend on this.choice. We suppose now that X is not an eigenvector 
of any of the O; so that the code has M distinct words. It would be 
desirable to be able to choose an X of this sort to either minimize P, 
or to maximize d, the nearest neighbor distance. We have not seen 
how to solve either of these problems in general. A few words about 
them are in order. 

Consider first the problem of sacs X to maximize d. The squared 
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distance between X and X, is 
d’(X, X,) = |X — X; |? = 2 — 2X-0O;X 


a monotone decreasing function of the quadratic form X-O;X in the 
components of X. This form is the cosine of the angle between X and X; . 
Solution of the maximum nearest neighbor distance is equivalent to 
finding 

a = min max X-0;X (9) 

x t 

where the maximization over the matrices of G must omit the identity 
matrix. The quantity a is an invariant of the representation (is the 
same for every equivalent representation) and should ultimately be 
expressible in terms of properties of the group. The vector X which 
minimizes (9) is not unique: any word in the code generated by X would 
serve as well. 

Given G, we define two points X and Y on the unit sphere to be equiv- 
alent if one can be obtained from the other by an operation of SG. The 
surface of the sphere is thus divided into equivalence sets. A connected 
region on the sphere such that no two points in its interior are equiv- 
alent and such that every point on the sphere is equivalent to some 
point in the region will be called a fundamental region of G. The maxi- 
mum likelihood regions, ®; , associated with any {M, n} generated 
by § intersect the unit sphere in fundamental regions. These inter- 
sections are very special fundamental regions: they are convex and 
bounded by hyperplanes. . 

In attempting to minimize P, or maximize d it clearly suffices to 
consider initial vectors X restricted to some fundamental region. It 
is natural then to ask what fundamental regions are possible for a 
given G. 

The situation is complicated. For some groups, the fundamental 
region is completely determined (up to equivalence under the group 
operations, of course): for other groups only certain features of its 
boundaries are determined, or no points at all may be determined. 

For example, in the plane consider the group G, generated by the 
three matrices corresponding to reflections in three lines through the 
origin that make angles of 60° with each other. This group is of order 6 
and is a subgroup of the symmetry group of a regular hexagon having 
the given lines as diagonals. The fundamental region of this group is 
completely determined. It is a 60° arc of the unit circle with end points 
on two of the given lines. Any group code {6, 2} generated by S, has 
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this fundamental region for the intersection of one of its maximum like- 
lihood regions ®; with the circle. Choice of X serves only to position 
the initial vector within the maximum likelihood region. (When X is 
chosen to lie on one of the reflection lines, a {3, 2} results and the maxi- 
mum likelihood region changes discontinuously to the union of two 
adjacent regions of the sort just discussed.) 

On the other hand, consider the group S, of rotational symmetries 
of the regular hexagon. G, , of order 6, consists of a 2 X 2 matrix rep- 
resenting a rotation of 60° in the plane along with the distinct powers 
of this matrix. Any 60° arc of the unit circle is a fundamental region for 
this group. Codes {6, 2} generated by SG, are equivalent for all choices 
of the initial vector X. 

An example illustrating a partly determined fundamental region is 
obtained by considering the pure rotational symmetries of a cube 
in three dimensions. We imagine the cube centered at the origin and 
inscribed in a unit sphere. We speak in terms of the operations on the 
cube rather than in terms of the 3 X 3 matrices which describe these 
operations. G;, a group of order 24, consists of rotations of the cube 
by 120° around the body diagonals, of rotations by 90° about axes 
through the origin and centers of faces and of rotations of 180° about 
axes through the midpoints of edges and the origin. One axis of each 
kind is shown on Fig. 1. In discussing the fundamental region of §; 
and codes generated by G;, it is convenient to speak of points on the 
cube, rather than on the circumscribed unit sphere. It is to be understood 


90° Ae A 180° 
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Fig. 1 — Example of partly-determined fundamental region. 
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then that when a point on the cube is mentioned it is really the cor- 
responding point on the sphere obtained by projecting along a radius 
that is under discussion. 

The vertices of the cube, the centers of faces and the midpoints of 
edges must all lie in boundaries of fundamental regions, for these points 
are on axes of rotation of G;. For example, a point distance e from a 
vertex of the cube has two nearby equivalent points forming an equi- 
lateral triangle with the vertex at the center of the triangle. These 
three points cannot all lie in the interior of one fundamental region. 
The cube vertex therefore cannot be an interior point of a fundamental 
region. In fact at least 3 fundamental regions must meet at each vertex, 
at least 4 at each face center and at least two at each edge midpoint. 
Cube vertices and face centers must therefore be vertices of fundamental 
regions. Now all vertices of the cube are equivalent under GS; as are all 
face centers and all edge midpoints; no two of these three types of 
points are equivalent. A fundamental region of G; must therefore con- 
tain at least one cube vertex and one face center among its vertices and 
at least one cube edge midpoint along its boundary. 

Two distinct types of fundamental regions for G; bounded by hyper- 
planes (great circles on the sphere) are shown in Fig. 1. Region AE FG 
is bounded by four hyperplanes. Edge midpoints are vertices of this 
type or region. Four fundamental regions surround each face center 
and each edge midpoint: three surround each cube vertex. Region 
ABCD is bounded by only three hyperplanes. Edge midpoints are 
no longer vertices of the fundamental region. Eight regions meet at 
each face center. The fundamental region ABCD corresponds to the 
maximum likelihood region of a group code having an initial vector 
(and hence all vectors) pass through a cube edge: region AEFG results 
when the initial vector passes through a face diagonal. All other positions 
of the initial vector give maximum likelihood regions that are funda- 
mental regions bounded by four hyperplanes but not congruent to 
AEFG. 

Gs is the irreducible representation of the symmetric group on four 
letters derived from the Young tableau”’ associated with the partition 
(2, 1, 1). The irreducible representation belonging to the partition (3, 1) 
is also three dimensional. It is equivalent to the group of symmetries 
of the regular tetrahedron and can be generated by reflections in planes 
through the centroid of the tetrahedron and its edges. The fundamental 
region here is completely determined. It is bounded by three of these 
generating reflection planes. Maximization of nearest neighbor distance 
for a {24, 3} generated by this group can be easily accomplished by 
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choosing the initial vector equidistant from the three bounding planes 
of the fundamental region. 

More generally, Coxeter” has shown that if a real-irreducible finite 
group of nm X n orthogonal matrices is generated by reflections, the 
fundamental region is completely determined, and in fact the region is 
bounded by 7 hyperplanes. Indeed, Coxeter has enumerated all possible 
groups of this sort. In dimensions n greater than 8, there are only three 
such groups, called by him A,,, B, , and C, of order (n + 1)!, 2"-*n!, and 
2°n!, respectively. These groups generate permutation modulation 
codes’—A, generates Variant I codes, B, generates Variant II codes 
with uw, = 0, and C,, generates Variant II codes with », ~ 0. The various 
permutation modulation codes are obtained by choosing the initial 
vector to lie in boundaries of various dimensionality of the fundamental 
regions of these groups. 

Returning to the general case (when G is not generated by reflections), 
the real eigenvectors of the O; with eigenvalue unity serve to determine 
landmarks of the fundamental region. Such an eigenvector must lie in 
the boundary of the region. If O; has / such eigenvectors, their span is 
an l-dimensional boundary of the fundamental region. The situation 
has been studied by Robinson™ in some detail, but no simple method 
of classifying the possible regions is available. 


V. THE DIRECT SUM 


Since any group code is equivalent to the weighted direct sum of codes 
generated by real-irreducible representations of a group, it is natural 
to investigate the relationship between interword distances in the sum 
code and the corresponding distances in the summand codes. 

Let G = A,, A., --:, A, be a finite group of order g with A, the 
identity. Let D'(A) and D*(A) be two real-irreducible representations 
of G by real orthogonal matrices of dimensions /, and /, respectively. 
Let X; = D’(A,)X, andY; = D*(A,)Y, 7 = 1, 2, --- , g be group codes 
generated by D* and D’. (Neither code need have g distinct vectors.) 
The direct sum code with weights \, and ), has vectors 


Z; = AX; P AY; = 1, 2: aie tg (10) 
pe ee ae 0O<rA,rAx2< l 


of 1 = 1, + l, components. We seek to choose the weights so that the 
nearest neighbor distance, d, for the sum code Z is a maximum. 

Let a; = d’(X;, X,) and 6; = d’(Y;, Y:) be the squared distance 
from the code word generated by A; to the initial vector in the two codes, 
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z= 1,2, --- , 9, respectively. For the sum code we have 
a (Z; » Zi) = Nia; — N38. 


since the subspaces containing the X code and the Y code are orthogonal. 
The desired maximum nearest neighbor distance is thus 

d’ = max min [(1 — dja; + 8] (11) 

OS\S1 141 

where we have set \ = 3. The situation is illustrated in Fig. 2. Here we 
have taken a. S a3 S +++: S a, which we can do without loss of general- 
ity since this is merely a matter of giving names to the group elements. 
The bracketed expression on the right of equation (11) is plotted as the 
line segment with ordinate a; at \ = 0 and ordinate 6; at \ = 1. We 
seek the highest point on the bottom boundary of this collection of 
lines, point P in Fig. 2. 

Now any of the vectors Y; ,7 = 1, 2, --- , g, not just Y, , would serve 
to generate the Y code. We can seek a further maximization of the nearest 
neighbor distance (11) for the Z code by choice of the vector from the 
Y code to be called Y, . Stated otherwise, for the initial vector of the Z 
code we choose a particular vector X, from the X code and to this we 
can add (directly) any of the vectors of the Y code. Now replacing Y, 
by Y, merely amounts to permuting the subscripts on the @; of Fig. 2. 
The subscript z is replaced by k where A,;A; = A;,. To combine the 
two codes to get the largest nearest neighbor distance, we must further 





0 1 X 


Tig. 2— Maximum nearest neighbor distance. 
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maximize (11) over permutations of the 6’s corresponding to left 
translations of the group. 

The maximization just considered was with respect to the manner of 
combining the two summand codes. There remains the matter of choos- 
ing X, and Y, to further increase (11). At first it might be thought 
that these vectors should be chosen to maximize the nearest neighbor 
distance in each of the summand codes. That this is not necessarily 30 
can be seen from Fig. 2. Choosing X, to increase the nearest neighbor 
distance in the code generated by D* would cause a, to increase. The 
line connecting a, and 6, on the figure would move up. However, this 
change of X, might cause a, to decrease by a larger amount so that 
point P on the figure moves downward. The situation is complicated. 

The relationship between the maximum likelihood region for the 
sum code and the corresponding regions for the constituent codes is 
even more complicated in general. Let ® be the region belonging to 
Z, of equation (10) and let ®’ and @’ be corresponding regions for X, and 
Y, in the summand codes. We write Z = \,X + d.Y fora general point in 
the space of the direct sum representation where X and Y lie in the 
respective invariant subspaces of the summand codes. A point will lie 
in ® then if |Z — Z, | < |Z — Z; | for 7 = 2, 3, --- , g, or what is the 
same, if 


Nid? (X, X,) ae ha’ (Y, Y;) = hid (K, X;) Ss a (Y, Y;) 


fori = 2,3, ---,g. Thusif X e@' andY e ® thenZ e ®, but the con- 
verse is not necessarily so in general. 

A special case in which the converse holds is the following. It may 
happen that both the X code and the Y code have fewer than g distinct 
vectors. In the direct sum code (10) it may happen that each distinct 
vector of the Y code is paired at least once with each distinct vector 
of the X code. (G must be homeomorphic to the direct product of two 
groups.) In this case ® is the cartesian product of the two regions ®* 
and ®*. The probability of no error for the sum code is given by Q, = 
Q3(\1)Q2(A2) where the factors are the probabilities of no error for the 
separate scaled summand codes. The information rate (3) for the sum 
code in this case is the weighted sum of the rates for the constituent 
codes 


hy 


ee bs 
R=FR, +78. 


We are better off using the code with the larger rate uncombined. 


GROUP CODES 589 
VI. THE CONFIGURATION MATRIX 


Let X,, X,, ---, Xa, be a collection of unit vectors spanning §, . 
The configuration matrix of this code is the M by M matrix p whose 
elements p;; = X;-X; are the cosines of the angles between the words. 
Equivalent codes have identical configuration matrices except for a 
possible relabeling of rows and columns. This configuration matrix is 
real, symmetric, non-negative definite and of rank n. The diagonal 
elements are unity and the off-diagonal elements are of magnitude no 
greater than unity. 

Conversely, we have the following 


Lemma: Every real, symmetric, M by M non-negative definite matrix of 
rank n with diagonal elements unity and off-diagonal elements of magnitude 
less than unity is the configuration matrix of a code of M unit vectors that 
span &, . 


The proof of this lemma follows readily from the fact that a real sym- 
metric M by M matrix p can be diagonalized by an orthogonal matrix O, 
that is, Op0~* = A where, since p is non-negative definite and of rank n, 
A has n positive diagonal elements and all other elements are zero. 
Without loss of generality we can take the first n diagonal elements of 
A, say \;;=); , 7=1, 2,---, n to be the positive ones. From p=O7*AO, it 
follows that 


Pii = by Ona Oui V Xu = X;-X; 


where X; is a vector of m components, the uth component being VO ui ‘ 
a = 1, 2, ---, M. We have now exhibited M unit n-vectors whose 
configuration matrix is the given matrix p. We need now only show that 
they span $, . But we have written p = XX where X is the matrix of M 
columns and 7 rows whose 7th column is X; . The tilde denotes trans- 
pose. Since the rank of a product of matrices is not greater than the 
smaller of the ranks of the factors, it follows that X must be of rank n, 
for if it were of rank less than n, so also would be p contrary to hypothe- 
sis. The X,; therefore span §, . 

For group codes, the rows of the configuration matrix are all permuta- 
tions of the first row of the matrix as can be seen from Proposition 1. 
Indeed the structure of this matrix is closely related to the multiplica- 
tion table of the group generating the code. Let the code vector X; = 
D(A;)X, 7% = 1, 2, --- , M be generated by an orthogonal representation 
D of a group G with elements A, , A2, --- , Aa, . Here A, is the identity 
and the code need not have M distinct vectors. Denote by @(A;) the 
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angle between X; and X, . Then 6(A;’) = 0(A,),7 = 1, 2, --- , M and 
the configuration matrix of the code is found to be given by 


pi; = cos 0(A;"A;) 


1,j = 1,2,---, M. If p,; < 1 for 7 > 1, then the code has JM distinct 
vectors: if 1 = pj, = py, = °° = prz, With 1 = 7, < jg ee- <j, 
and these are the only elements of value unity in the first row, then 
the code has M/h distinct vectors. 

Conversely, we have 


Theorem 1: Let x(A;) be a real-valued function defined on the elements 
Ai,As, ++: ,An, of a group G of order M. Let x(A,) = 1, where A, ts the 
identity of the group, and let x(A;) = x«(A;"),j = 1, 2, --- , M. If the 
M by M matrix p with elements p;; = x(Az'A;) is non-negative definite 
and of rank n, then there exists a group code {M’ n} generated by an n- 
dimensional orthogonal representation of G that has configuration matrix p. 
Here M’ = M/h where h is the number of different values of j for which 
x(A;) = 1. 


Proof: The proof follows easily from the lemma. We can find M unit 
vectors X; (not necessarily distinct) that span §, such that p;; = X;-X; . 
Without loss of generality we can suppose that X,, X., --- , X, are 
linearly independent. Now an n by n real matrix is determined by speci- 
fying its effect on n vectors that span its carrier space. For each np = 
1, 2, --- , M we determine the n by n matrix D(A,) by specifying its 
effect on X,,---, X,, namely that D(A,)X; = Xigjy,7 = 1,2,---,n 
where A,A;=Ai(,,). Now X;-X;=p;;=2(Ajz"A,;) =a2(A;"A7'4,4;) = 
UAT wAriw) = Xrw*Xi) , 80 that D(A,) preserves the angles 
between n vectors spanning its carrier space. It is easy then to show 
that D(A,) preserves the angle between any two vectors and hence is an 
orthogonal matrix. For 7 > n, 


X;= >»; ajnXn 
h=1 


for some set of a’s. Using this representation and the orthogonality of 
D(A,), it is now easy to show that D(A,)X;= Xi,;,,. for 7=1, 2, --- , M. 
The fact that D is a representation then follows readily. 

Theorem 1 permits an interesting reformulation of the problem of 
finding an {M, n} of largest nearest neighbor distance generated by a 
representation of §. Form the modified multiplication table of S,—an 
M by M array of group elements with A;*A; in the 7th row and jth 
column. From this table we construct a symmetric M by M matrix p 
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by replacing the group identity A, , say, by unity, by replacing both A, 
and Ay’ by the variable z,, A, and Az’ by x,, and so on. If G has 
exactly m self-reciprocal elements, there will be K = m — 1+ (g — m)/2 
variables in p. The condition that p be non-negative definite and of rank 
not greater than n obtained by conditioning certain minors of p gives 
rise to polynomial constraints in the variables z,, --- , xx. To find the 
code of largest nearest neighbor distance, we must minimize max; 2; 
subject to these constraints. 

We notice in closing this section that the configuration matrix of an 
{M, n} generated by a group G of order MZ commutes with all the ma- 
trices of the regular representation of G (see the appendix). Using Schur’s 
lemma, one can then arrive at a canonical representation for configura- 
tion matrices that involves the irreducible representations of S. But 
we do not pursue this topic further here. 


VII. SOME THEOREMS ON DISTANCES 


We now adopt the notation of the appendix. Let G be an abstract 
group of order g with elements H, A, B, --- where E is the identity. 
Let D(Z), D(A), --- bea real-irreducible representation of G by n X n 
(real) orthogonal matrices. From an initial unit vector K = X, the 
representation generates a code by Xz, = D(R)X,, RF runs through S. 
We denote the squared distance from X, to Xs by d’(KXz, Xs). We 
have then 


e-€ , x) — yy —™ 2 De D(A) 50,2; 
t,j=1 
— a (ax , Xr) (12) 
for every R and A ¢G. Here x,, 22, --* , x, are the components of X. 
For codes generated from real-irreducible representations in this 
manner, a number of interesting distance sums are independent of the 
choice of the initial vector X. 


Theorem 2: Let D(R) be the matrices of a real-irreducible orthogonal repre- 
sentation of a group of order g. Let Xp = D(R)X. If DUR) ts not the trivial 
one-dimensional representation D(R) = 1, then 


2 d’(Xp , X) = 2g 
ReG 
independent of the unit vector X. 


This is really a special case of the more general 
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Theorem 3: For any code generated from the initial unit vector X by a 
real-irreducible orthogonal representation D of G, 


» d’(Xpm ? X) = 2g(1 > Lm) 
where 


gee) 


gn Reg 


as a constant independent of X. Here x(R) = Tr D(R) ts the character of R 
in the representation. 


Proof: Consider the matrix 
Az » D(R”) = p> D(R)D(R) --- D(R) 


where there are m factors in the summand. Since the representation is 
by orthogonal matrices, D(R) = D™*(R) = D(R™*) where the tilde 
denotes transpose. Thus 


A = > D(R”)--- DR) =A 
Re§ 


since as R runs through § so does R™*. The matrix A is thus symmetric. 

We next show that A commutes with all the matrices D(R). By a 
theorem quoted in the appendix we can then conclude that A = al 
where I is the unit matrix. To see that A commutes with D(R), consider 


AD(R) = 2, D(S)""D(S)D®) = p> D(S)""D(SR). 


Now set SR = T's0 that S = TR™*. Then 
AD(R) = >> D(TR™)""D(T) 
TG 


>> D(TR™)D(TR™) --- D(TR™)D(T) 
> DL)D(R"T)D(RUT) --» DRT) 
T eG 


> D(RU)D(U)"* = DR) 2, DU)” = D(R)A 
UeG Us 


where we have used the substitution U = R™’T. 
From equation (12) we have 


4,7=1 


> d’(Xpn ,X) = 2g — 2 >; Uh; oy D(R");; 
ReG ; ReG 


2g — 2 De ve;Ai; = 2(g — a) 


t,7=1 
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by the diagonal property of A just established. To find a consider the 
trace of A. We have 


Tr A = Tral =an = >> Tr D(R”) = > x(R”). 
ReG 


ReG 


The theorem then follows. 

To establish Theorem 2, notice that for the trivial representation 
D'(R) = 1, we have x’(R) = 1. For the character x(R) of any other 
nonequivalent real-irreducible representation we then have 


» x'(R)x(R) = p> x(R) = 0 


by the orthogonality relations (appendix). Using this fact and setting 
m = 1 in Theorem 8 yields Theorem 2. 


Theorem 4: Let © be a class of n, elements of G with character x(C). For 
any code generated from the initial unit vector X by a real-irreducible 
orthogonal representation D of S, 


E Ke ,¥) = ani — 20) a3) 


‘Ree 

| independent of the unit vector X. 

Proof: 
>. d’(Ke , X) = 2n. — 2 Do az; D5 DR); (14) 
Ree Re 


Now consider the matrix 


B = 3) DR) =" Y DSRS) = ™ DY D(8)D(R)D(S") 
Ree J seg J seg 

where D“(R) is an irreducible (over the complex field) representation of 

dimension m of S. Now B commutes with all the matrices of D* since 


BD°(T) = 7 > D°(S)D°(R)D*(S"P) 
a >> DTU) D*(R)D*(U) = D*(T)B 
UeG 


where we have set S-'7' = U. By Schur’s lemma, B = kI where I is 
the m by m unit matrix. Taking traces we have 


TrB =Tr >> D*(R) = 2.x°(C) = Trh&l = km 


Ree 
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so that 


B= DDR) =" x"(OL. (15) 
Re€ m 

If now the real-irreducible orthogonal representation D is also irre- 
ducible, by equation (15) the inner sum in (14) is (n,/n)x(C)6;; and the 
theorem (13) follows at once. 

Suppose now that D is not irreducible. Then (see appendix) D is 
equivalent to an orthogonal representation of the form 


UR) V"(R) 
—V*(R) U*(R) 


where D°(R) = U*(R) + iV“°(R) is an irreducible representation by 
unitary matrices and U* and V°* are real and of dimension m where 
n = 2m. We can suppose the D of equation (14) to be of the form (16). 
Now let 


DR) = (16) 








B= >) DR) 
Ree 
and set 


-7) 
y 1 ae 


V2 lr 7 
where as before I is the m by m unit matrix. One then finds by direct 
computation that 











- D‘(R) 0 

U"BU = > 
Ree 0 D*(R)* 
me |x(C}L 0 | Ly 
™L 0  x"@*L 








where the middle equality follows from equation (15). Now let x*(@) = 
pw + w with yp and v real. Direct computation gives 


wl vl) 
—vl pl 


(17) 





B = UHU? = | 
nm 


The right of (14) is 
2n, — 2 25 B, 2,2; 
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and using (17) this becomes 


2n. — iL ud t= an,(1 a H), 


mM 1 





From equation (16), however, x(C) = 2u, so that (13) then follows 
and the theorem is proved in all cases. 

Since every group code can be thought of as the direct sum of codes 
generated by real-irreducible representations, and since squared distance 
in the sum code is the sum of squared distances in the separate codes, 
Theorems 2, 3, and 4 have ready analogues for all group codes. For 
example, if a group code does not contain the identity representation, 
then 

Da (Xe , X) = 2g. 

ReG 
Theorems 3 and 4 hold for group codes in general when x(R) is replaced 
by the weighted sum >) \*x‘(R) of the characters of the constituent 
real-irreducible codes. 

Another theorem of interest concerning codes generated from any 
group of orthogonal matrices arises from the fact (12) that d’(X,, X) = 
d’(Xra, Xx). Let there be a point of the code distant d from the point 
X,. Starting from Xz, we imagine moving from word to word of the 
code restricting our moves so that from any word we can move only to a 
word distant d away. We shall call the collection of words that can be 
reached from X, in this manner “‘ad chain starting from X,’’. Xz is to be 
included in this chain. 


Theorem 5: Let the words of a d chain starting from Xzg be Xz, Xa,, 
Xu,, °°: Xa,- Then the group elements i, Ay, Az, °+: , A, forma 
subgroup 3C of G. The group elements whose corresponding words are dis- 
tant d from Xz form a set of generators for 5C. If 5C is a proper subgroup of G, 
then from any word corresponding to a group element not in 3C, a new d 
chain may be formed and the group elements corresponding to the points of 
this new d chain will form a coset of 3. 


Proof: Suppose all the points distant d from Xz are X4, , Xa,,°°-, Xa, - 
Let us construct a table of group elements in the following manner. 
The first row is H, Ay, An, ---,A,. The K + Ist row of the table is 
formed from the preceding K rows as follows. We examine the elements 
of the table in order, reading from left to right in the first row, then from 
left to right in the second row, and so on. Let R be the first element 
arrived at in reading the first K rows that does not appear in the first 


596 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1968 


column, rows 1, 2, --- , K. The the K + Ist row is to be 
R, RA, 5, Ag, eee , RA, . 
The table thus appears 


EK, A, ) A, ) Ass 
A, ? : ) A, Az ) A, An 


Ay ) A.A, 5) Az ) AoA 


Aes, Males “dix 2 
B, BA, , BAS BA, 


R, RA, , RAz , RA, 


When j rows have been written and every element in these 7 rows has 
appeared once in the first column the process is stopped and the table is 
considered complete. The table can have at most g rows. Now from 
d?(Xz, Xz) = d’(Xsz, Xs), it follows that the words represented by 
the elements in the 2nd, 3rd, --- , m + 1st columns of the Kth row are 
all distant d from the word represented by the element in the first column 
of the Kth row. Furthermore, these m words are all the words of the 
code that are distant d from the word represented by the element in the 
first column of the Kth row. Thus the elements of the first column of 
the table give the points of the d chain starting from Xz. That these 
elements of the first column form a group 3C and that A, , Az, -::,An 
are generators of 3C is clear from the method of constructing the table, 
for we have formed all possible distinct products of the A’s and listed 
the distinct elements thus obtained in the first column. Let 3C be a 
proper subgroup of G and let S be an element of G not in &. If we multiply 
every element in the above table by S, we obtain a new table giving 
all the points that can be reached from point Xs by steps of distance d. 
The first column of this table lists the points of the d chain starting 
from Xs and the corresponding elements are just the coset S3c of 3c. 


VIII. BINARY GROUP CODES 


The group codes (n, k) for the binary channel introduced in Ref. 21 
are group codes in the present sense. Each word of an (n, k) code is an 
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n-place binary sequence. Replace each zero by 1 and replace each 1 by 
—1 in each word. Then write each word (a sequence of +1’s) as a 
diagonal n X n matrix. This collection of 2° n X n orthogonal matrices 
forms an Abelian group @, that is isomorphic to the k-fold direct product 
of the simple two element Abelian group. The matrices generate the 
code by operating on the n-vector (1, 1, 1, --- , 1). The real-irreducible 
representations of this group are all one dimensional. There are 2" of 
them. The representation by n X n matrices just considered is already 
exhibited in reduced form as the direct sum of n of these real-irreducible 
representations. 


IX. CONCLUDING REMARKS 


The foregoing paragraphs outline some of the theory of group codes 
for the Gaussian channel. The development of this subject is clearly 
incomplete: we have raised more questions than we have answered. 
Perhaps the outstanding problem is that of finding a tractable method 
of choosing the initial vector to maximize the nearest neighbor distance. 

There is a great abundance of groups of arbitrarily large order that 
can be examined from the point of generating group codes. The sym- 
metric group and the hyperoctahedral group appear most promising 
for initial investigation since their structure and irreducible representa- 
tions (which are all real) are comparatively well understood. 


APPENDIX 


Review of Group Representation Theory” 


Let G be a finite group of order g with elements H, A, B, --- . The 
letters R and S will be used for the general element of and EF will denote 
the identity of §. As R runs through G, the distinct elements of the set 
RAR ™ are said to form a class of §. The elements A and B are said to 
belong to the same class of ¢ if there exists an S such that A = SBS™. 
G can be divided uniquely into a union of classes, no two classes con- 
taining a common element. The number of elements in a class of G is a 
divisor of g. 

If 3C is a subgroup of G and if 3C is of order h, then h is a divisor of 
g and the number g/h is called the index of 3¢ under S. If, for every 
R in 3, all elements of G in the same class as Ff are also contained in 3, 
then 3C is said to be a self-conjugate subgroup of S. A subgroup JC of S 
is sald to be proper if h < g. 

The matrices in what follows are assumed to have elements in the field 
of complex numbers. 


598 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1968 


If to every element F of a finite group § there corresponds an n by n 
nonsingular matrix D(F) and if D(R)D(S) = D(RS), the collection of 
matrices A = {D(R), R runs through G} is said to form an n-dimensional 
representation of G. The matrices of A form a group under matrix 
multiplication. If the correspondence between the matrices of A and the 
elements of G is one-to-one, A is said to be a faithful representation of S. 
If for some R ~ S, D(R) = D(S), A is said to be an unfaithful repre- 
sentation of S. The matrix D(E) is always the n by n unit matrix. If a 
representation is unfaithful, the elements represented by D(E) form a 
self-conjugate subgroup of S, say of order h, and to each matrix of A 
correspond exactly h elements of §. A contains g/h distinct matrices. 
If D(Z), D(A), --- is an nm dimensional representation of §, so is 
MD(E)M~, MD(A)M~, --- where M is any nonsingular n by n 
matrix. The two representations A and MAM™ are called equivalent. 
Every representation of a finite group is equivalent to a representation 
by unitary matrices. Henceforth we shall be concerned only with such 
unitary representations. 

A finite collection of n by n matrices 0, , O2, --- , Ox is said to be 
reducible if there exists an n by n unitary matrix U such that for 7 = 
1,2, --- , K we have 


D 


B; 


Ay 
C 








U0;U~ = | 


where A; is an l by J matrix, B; is an n-l by n-l matrix,0 <1 <n,C 
is an n-l by | matrix all of whose elements are zero, and D is an l by n-l 
matrix all of whose elements are zero. It is assumed that / is independent 
of z. A collection of matrices that is not reducible is said to be irreducible. 

Every finite group has exactly as many nonequivalent irreducible 
representations as it has classes. If 1, ,1,, --- , J, are the dimensions of 
all the nonequivalent irreducible representations of S, of order g, then 


2 Sa, 
1 


If D*(R),, is the element in the wth row and vth column of the matrix 
representing FR in the /,-dimensional irreducible representation, a, of S, 
then 


Ds D*(R),,D*(R)E,- = Bunt b9/la My, b's Vy y= Ls 2, ea ly.3 
ReG 


Here * means complex conjugate and 6 is the usual Kronecker symbol. 
If the matrices D*(R) form an 1, dimensional irreducible representation 
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of G not equivalent to the representation a, then 


Dd, DR)» DBE, = 


ReG 
py Ly ep SOR dean fe gl Vy 2h. etry dees 


If D(R) is the n by n matrix representing F# in the representation A, 
the trace of D(R), namely 


xR) = Y Dan, 


is called the character of R in the representation A. If R and S are in 
the same class of G, then x(R) = x(S), for any representation of G. 
The characters of the irreducible representations a and 6 of G satisfy 
the orthogonality conditions 


2d x"(R)x(R* = 9 Sa - 


Here 6,., is unity if a and 8 are equivalent representations and is zero 
otherwise. 

Let A be any representation of S with character x(R). Let the char- 
acters of the irreducible representations of G be x’(R), 7 = 1,2, --+,¢ 
where c is the number of nonequivalent irreducible representations of 
GS (= number of classes of S$). Then x(R) may be written uniquely in the 
form 


xR) = DYax(R), allRing, 


where the a; are nonnegative integers independent of R. In fact, 


a, = 7 xR)", 
ReG 


The representation A is said to contain the irreducible representation 
j a; times and there exists a unitary matrix U independent of R such 
that 


D*(R) ee 0 
UD(R)U™* = 0 Di(R) --- 0 all Rk ng 


0 0 DR) 


600 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1968 


where D*(R), D'(R), D*(R) etc., are matrices of the ith, jth, kth, etc., 
irreducible representation of G, and 0 stands for the appropriate matrix 
with all elements zero. The jth irreducible representation will occur 
exactly a; times among D‘(R), D’(R), D*(R) and so on. 

Every group § possesses a faithful representation, called the regular 
representation I, that consists of g by g permutation matrices. The rows 
and columns of these matrices can be labelled by the elements of G. The 
entry in row & and column S of the matrix representing T is unity if 
Rk = TS and is zero otherwise. The regular representation is reducible: 
it contains the irreducible representation D*“ exactly /, times, a = 
BD, ekg 

Let D’(R) and D’’(R), R runs through G, be irreducible representations 
of S of dimension d@’ and d” respectively. Let the matrix H satisfy 
D'(R)H = HD"(R) for all R in . Then either H is the zero matrix or 
H is square and nonsingular so that d’ = d’’, and the two representations 
are equivalent. A matrix that commutes with all the matrices of an 
irreducible representation of G is a multiple of the unit matrix. These 
statements are known as Schur’s lemma. 

Much of the foregoing remains valid with minor modifications when 
the number field in question is the real rather than the complex numbers. 
One easily finds that every real representation of a finite group is 
equivalent (over the reals) to a representation by orthogonal matrices. 
The only real symmetric matrix that commutes with all the matrices 
of a real-irreducible representation is a multiple of the unit matrix. 
If D*(R) and D*(R) are nonequivalent real-irreducible representations 
by real orthogonal matrices, respectively of dimension J, and [,, then 


2s D*(B)wD°(R) ary a 0, 
ReG 
B,p = 1,2, 0%9 ole ye dy 25-8 ley 


» LD*) yD" (BR) wry T D*(R) yy D*(B) ys) = 2 Duy’ Sy G/ ba ’ 
ReG 


BV ww SH 1,2, be. 
For the characters one has 


>, x*(R)x(R) = 0 
ReG 
if the representations a and 6 are not equivalent, while 


a [x*(R)x*(R) + x*(R’)] = 29. 
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Every real-irreducible representation that is reducible (over the complex 
numbers) is equivalent to the direct sum of an irreducible representation 
and its complex conjugate. If the irreducible unitary representation 
D(R) = U(R) + iV (Rk), with U and V real, is not equivalent to a real 
orthogonal representation, then 


U(kt) | V(R) 
(ae ie) (18) 


is a real-irreducible representation by real orthogonal matrices. 
For an irreducible representation D(R) with character x(R), the sum 


n=t SR) 
g ReG 
can have only one of the three different values 0, +1. If h = 1, D(R) is 
equivalent to a representation by real orthogonal matrices. Ifh = —1, 


the representation D is equivalent to its complex conjugate, but is not 
equivalent to a real representation. A real-irreducible representation 
can be made from each irreducible representation D having h = —1 by 
forming real matrices of the form (18), where U and V are the real and 
imaginary parts of D. Finally, if h = 0, D is not equivalent to its com- 
plex conjugate and is not equivalent to a real representation. Non- 
equivalent irreducible representations for which h = 0 occur then in 
complex conjugate pairs. Each such pair gives rise to a single real- 
irreducible representation through the recipe (18). Thus, finally, if h 
has the value 0 for exactly 2p of the c nonequivalent irreducible repre- 
sentations of G, then G has exactly c — p nonequivalent real-irreducible 
representations. 


REFERENCES 


= 


Kotel’nikov, V. A., Thesis, Molotov Energy Institute, Moscow, 1947, trans- 
lated as The Theory of Optimal Noise Immunity, New York: McGraw-Hill 
Book Co., 1959. 
. Shannon, C. E., “A Mathematical Theory of Communication,” B.S.TJ., 27, 
Nos. 3 and 4 (July and October 1948), pp. 379-423; 623-656. 

. Shannon, C. E., “Communication in the Presence of Noise,” Proc. IRE, 387, 
(January 1949), pp. 10-21. 

. Rice, 8. O., “Communication in the Presence of Noise—Probability of Error 
for Two Encoding Schemes,” B.S.TJ., 29, No. 1 (January 1950), pp. 60-93. 

. Gilbert, E. N., “A Comparison of Signalling Alphabets,” B.S.TJ., 31, No. 3 
(May 1952), pp. 504-522. 

. Shannon, C. E, “Probability of Error for Optimal Codes in a Gaussian 
Channel, "BST J., 88, No. 3 (May 1959), pp. 611-656. 

. Wolfowitz, J., Coding Theorems of Information Theory, Berlin: Springer- 
Verlag, 1961. 

. Stutt, C. A., “Information Rate in a Continuous Channel for Regular- 

Simplex Codes,” IRE Trans. JT-6 (December 1960), pp. 516-522. 


ont Om oO Fe WW WW 


602 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1968 


9. 


10. 
i; 
12. 
13. 
14. 


15. 
16. 


17. 


18. 
19. 
20. 
21. 
22. 
23. 
24, 


25. 


Balakrishnan, A. V., “A Contribution to the Sphere-Packing Problem of 
Communication Theory,” J. Math. Anal. Applications, 3 (December 1961), 
pp. 485-506. “Signal Selection Theory for Space Communication Chan- 
nels,” Chapter 1 in Advances in Communication Systems, ed. A. V. . 
Balakrishnan, New York: Academic Press, 1965. 

Slepian, D., “Bounds on Communication,” B.S.TJ., 42, No. 3 (May 1963), 
pp. 681-707. 

Gallager, R. G., “A Simple Derivation of the Coding Theorem and Some 
Applications,” IEEE Trans., /7-11 (January 1965), pp. 3-18. 

Wyner, A. D., “Capabilities of Bounded Discrepancy Decoding,” B.S.TJ., 
44, No. 6 (July-August 1965), pp. 1061-1122. ; ; 

Wozencraft, J. M. and Jacobs. I. M., Principles of Communication Engi- 
neering, New York: John Wiley & Sons, 1965. 

Weber, C. L., “On Optimal Signal Selection for M-ary Alphabets with Two 
Degrees of Freedom,” IEEE Trans., IT-11 (April 1965), pp. 299-300, and 
“New Solution to the Signal Design Problem for Coherent Channels,” 
IEEE Trans., [7-12 (April 1966), pp. 161-167. 

Slepian, D., “Permutation Modulation,” Proc. IEE 53 (March 1965), pp. 
228-236. 

Zetterberg, L. H., “A Class of Codes for Polyphase Signals on a Band- 
limited Gaussian Channel,” IEEE Trans., J7-11 (July 1965), pp. 385-395. 
“Detection of a Class of Coded and Phase-Modulated Signals,” JTEEE 
Trans., [7-12 (April 1966), pp. 153-161. 

Peterson, W. W. and Kasami, T., “Reliability Bounds for Polyphase Codes 
for the Gaussian Channel,” Scientific Report No. 3 (July 1965), Dept. of 
Elec. Eng., U. of Hawaii. Abstract appears in IEEE Trans., [T-12, No. 2 
(April 1966), p. 277. 

Reed, I. S. and Scholtz, R. A., “N-Orthogonal Phase-Modulated Codes,” 
IEEE Trans., [T-12 (July 1966), pp. 388-395. 

Scholtz, R. A. and Weber, C. L., “Signal Design for Phase-Incoherent Com- 
munications,” IEEE Trans., 77-12 (October 1966), pp. 456-463. 

Landau, H. J. and Slepian, D., “On the Optimality of the Regular Simplex 
Code,” B.S.TJ., 46, No. 8 (October 1966), pp. 1247-1272. 

Slepian, D., “A Class of Binary Signaling Alphabets, B.S.T.J., 35, No. 1 
(January 1956), pp. 203-234. 

Dunn, James G., “Coding for Continuous Sources and Channels,” Thesis, 
Electrical Engineering Department, Columbia University, May 1965. 

Coxeter, H. S. M., Regular Polytopes, New York: The Macmillan Com- 
pany, 1963. 

Robinson, G. de B., “On the Fundamental Region of an Orthogonal Repre- 
sentation of a Finite Group,” Proc. London Math. Soc., 43, Sec. 2 (1937), 
pp. 289-801. 

Boerner, H., Representation of Groups, Amsterdam: North-Holland Publish- 
ing Co., 1963. 

Murnaghan, F. D., The Theory of Group Representations, Baltimore: 
Johns Hopkins Press, 1938. 

Weyl, H., The Classical Groups, Princeton: Princeton University Press, 
1946. 

Wigner, E. P., Group Theory, New York: Academic Press, 1959. 
a J.§., Applications of Finite Groups, New York: Academic Press, 


Contributors to This Issue 


Jack M. Hourzman, B.E.E., 1958, City College of New York, M5S., 
1960, University of California (Los Angeles) ; Ph.D., 1967, Polytechnic 
Institute of Brooklyn; Hughes Aircraft Company, 1958-1963; Bell 
Telephone Laboratories, 1963—. Mr. Holtzman has worked in various 
aspects of systems and control theory. Member, SIAM, Sigma Xi. 


Lupwik Kurz, B.E.E., 1951, and M.E.E., 1955, City College of New 
York, Eng. Se.D., 1961, New York University. Mr. Kurz is a professor 
of electrical engineering and project director in the Laboratory for 
Electroscience Research of the School of Engineering and Science at 
New York University. His major interests are in communication sys- 
tems optimization and application of the theory of stochastic processes 


to problems in electrical engineering. Member, Eta Kappa Nu, Sigma 
Xi. 


Rosert W. Lucky, B.S.E.E. 1957, M.S.E.E. 1959, and Ph.D. 1961, all 
from Purdue University; Bell Telephone Laboratories, 1961—. Mr. 
Lucky has been concerned with various analytical problems in the 
transmission of digital information over voice telephone facilities. 
He is Head of the Data Theory Department. Member, IEEE, Tau 
Beta Pi, Eta Kappa Nu, Sigma Xi. 


R. Maenanti, A.B., 1961, Columbia College; B.S.E.E., 1962, 
M.S.E.E., 1964, both from Columbia University; Bell Telephone Lab- 
boratories, 1962—. Mr. Magnani has worked on specifying system 
requirements for local, toll, and military electronic switching systems. 
He is a member of the Local Switching Studies Department and is 
engaged in statistical analysis and long-range forecasting. Member, 
Eta Kappa Nu. 


Bernarp T. Murpuy, B.Sc., 1958, Ph.D., 1959, both from the Uni- 
versity of Leeds, England; Bell Telephone Laboratories, 1963—. Mr. 
Murphy worked in the field of medical physics at the University of 
Leeds, on electron beam studies at Mullard Research Laboratories, 
and since 1959 has been engaged in work on semiconductor devices. 
At Bell Laboratories, he has worked on both the circuit and structural 
aspects of semiconductor integrated circuits. He supervises work on 


603 


604 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1968 


IMPATT diodes, high speed pulse generation, and new integrated 
circuit structures. 


Irwin W. Sanpperc, B.E.E., 1955, M.E.E., 1956, and D.E.E. 1958, 
Polytechnic Institute of Brooklyn; Bell Telephone Laboratories, 
1958—. Mr. Sandberg has been concerned with analysis of military 
systems, synthesis and analysis of active and time-varying networks, 
studies of properties of nonlinear systems, and some problems in 
communication theory and numerical analysis. He is head of the 
Systems Theory Research Department. Member, IKEE, Eta Kappa 
Nu, Sigma Xi, Tau Beta Pi. 


H. Suicuman, B.S., 1962, University of Michigan; M.S., 1963, Uni- 
versity of Illinois; Bell Telephone Laboratories, 1962—. Mr. Shich- 
man is a member of the Digital Device Integration Department where 
he has worked on characterization and design of logic circuits. He 
is working on computer analysis and synthesis of nonlinear networks. 
Memker, IEEE, Eta Kappa Nu, Tau Beta Pi, Phi Kappa Phi. 


Davin Stepian, 1941-48, University of Michigan; M.A., 1947. 
Ph.D., 1949, Harvard University; Bell Telephone Laboratories, 
1950—. He has been engaged in mathematical research in communi- 
cation theory and the theory of noise, as well as in a variety of as- 
pects of applied mathematics. Mr. Slepian has been mathematical 
consultant on a number of Bell Laboratories projects. During the aca- 
demic year 1958-59, he was Visiting Mackay Professor of Electrical 
Engineering at the University of California at Berkeley and during 
the Spring semester 1967 he was a Visiting Professor of Electrical 
Engineering at the University of Hawaii. Member, AAAS, American 
Math. Society, Institute of Math. Statistics, IEEE, SIAM. 


Nicuouas A. StraKHov, B.S.M.E., 1959, Massachusetts Institute of 
Technology; M.E.E., 1961, Ph.D., 1967, New York University; Bell 
Telephone Laboratories, 1959—. Mr. Strakhov has been designing 
and developing electronic test sets for transmission media mainte- 
nance. He supervises a group responsible for developing new trans- 
mission media. Member, Sigma Xi, Pi Tau Sigma, IEEE. 


